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Section  1 

Introduction 


Material  response  under  impact  loading  conditions  is  very  complex.  Most  often,  the 
impact  loading  induces  compressive  shock  under  one-dimensional  strain;  however,  later,  the 
reflected  waves  from  the  lateral  stress  free  boundaries  generate  intense  tensile  loading.  In  many 
applied  engineering  problems,  the  stress  and  strain  states  are  three-dimensional  and  the  strain 
rates  are  extremely  high.  The  material  deformation  and  failure  processes  are  greatly  influenced 
by  these  loading  conditions.  The  initially  intact  material  develops  microvoids  and/or  microcracks, 
which  lead  to  a  loss  of  stiffness.  The  coalescence  of  these  defects  ultimately  leads  to  complete 
failure  of  the  material  under  high  velocity  impact.  Mathematical  models  describing  these 
fundamental  damage  processes  are  important  tools  in  the  analysis  of  target  and  warhead  designs 
using  advanced,  shock-wave-based,  finite  element/difference  computer  codes.  With  the  advent 
of  supercomputers,  it  has  become  feasible  to  use  computationally  demanding,  but  physically 
based  constitutive  models  as  predictive  tools  in  impact  engineering  calculations. 

This  report  summarizes  the  development  of  advanced,  three-dimensional,  constitutive/ 
damage  models  for  metals  and  ceramics.  These  models  are  suitable  for  a  variety  of  dynamic 
loading  conditions.  The  model  development  philosophy  was  to  (1)  describe  most  of  the  physical 
processes  that  operate  in  the  material,  (2)  minimize  the  mathematical  complexity,  (3)  keep  the 
number  of  model  constants  to  a  minimum,  and  (4)  devise  an  efficient  numerical  scheme  to 
implement  into  any  general  purpose  shock-wave-propagation  based  computer  code. 

Under  this  program,  a  new  dynamic  model  (the  RDG  model)  to  describe  void  nucleation, 
growth,  and  coalescence  in  ductile  metals  was  developed  by  Rajendran,  Dietenberger,  and  Grove 
[1].  This  three-dimensional  constitutive  model,  based  on  a  pressure  dependent  yield  criterion  for 
compressible  plastic  flow,  is  strain  rate  and  loading  history  dependent.  A  detailed  description 
of  the  RDG  model  as  well  as  the  constants  determined  for  several  metals  are  reported  in  Section 
3.  The  model  was  incorporated  into  the  EPIC-2  code  and  several  example  problems  were 
successfully  simulated.  The  results  of  these  simulations  are  also  summarized  in  Section  3.  A 
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impact  response  of  such  altered  material.  In  hydrocode  calculations,  a  simple  linear  relationship 
between  the  strength  and  pressure  is  assumed  for  the  powdered  ceramic. 

Characterizing  the  impact  behavior  of  ceramic  materials  involved  the  following  impact 
experimental  techniques:  (1)  split  Hopkinson  bar  (SHB),  (2)  plate  impact,  (3)  plate  impact  on 
a  bar,  and  (4)  bar-on-bar  impact.  In  these  experiments,  strain  gauges,  stress  gauges,  and  high 
speed  cameras  were  used  to  measure  strain,  stress,  displacement,  and  deformed  shapes.  The 
results  of  these  experiments  were  used  to  calibrate  the  material  constants  in  the  model. 

In  addition  to  the  development  of  the  metal  and  ceramic  models,  this  program  involved 
extensive  material  characterization  under  dynamic  loading.  Much  of  this  information  was 
obtained  in  support  of  the  model  development.  The  high-strain-rate-test  data  obtained  under  this 
program  for  several  metals  and  ceramics  are  summarized  in  Section  2. 

Both  the  fragmentation-based  and  microphysical-based  models  have  been  incorporated 
into  the  EPIC-2  code.  Laboratory  impact  experiments  on  AD85  ceramics  were  simulated.  Using 
the  advanced  microphysical  model,  a  sample  calculation  of  a  long  rod  penetration  into  a 
confined  ceramic  plate  was  also  performed.  The  results  of  the  numerical  simulations  as  well  as 
a  complete  description  of  the  ceramic  model  are  presented  in  Section  4. 

The  summary  and  recommendations  are  discussed  in  Section  5.  Appendix  A  summarizes 
the  SHB  and  plate  impact  test  data.  The  Bodner-Partom  [13]  and  Johnson-Cook  [14]  strength 
models  are  presented  in  Appendix  B.  The  numerical  algorithm  and  the  solution  technique  to 
solve  the  RDG  model  are  discussed  in  Appendices  C  and  D,  respectively.  The  modified  TCK 
model  used  to  describe  the  impact  damage  in  ceramics  is  presented  in  Appendix  E. 
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Section  2 

High  Strain  Rate  Test  Results 


The  two  main  objectives  of  the  high  strain  rate  testing  are  (1)  to  determine  the  dynamic 
properties,  such  as  the  variation  of  strength  with  respect  to  strain  rate,  Hugoniot  elastic  limit 
(HEL),  spall  strength,  and  initial  yield  strength,  and  (2)  to  obtain  the  dynamic  deformation  and 
stress  characteristics  from  high  speed  photographs  of  dynamically  deforming  (metal)  and 
fracturing  (ceramics)  specimens  using  a  high  speed  camera,  and  stress-time  history  data  from 
stress  gauges. 

For  this  purpose,  several  impact  test  configurations  were  considered:  (1)  split  Hopkinson 
bar,  (2)  plate  impact,  (3)  flyer  plate  impacting  a  long  bar,  and  (4)  short  bar  impacting  a  long  bar. 
In  this  section,  no  attempt  is  made  to  describe  the  test  techniques  or  the  data  analysis;  however, 
brief  descriptions  of  split  Hopkinson  bar  and  plate  impact  tests  and  a  summary  of  data  are 
provided  in  Appendix  A. 

2.1  Split  Hopkinson  Bar  Data 

The  split  Hopkinson  bar  (SHB)  experiments  provide  stress-strain  data  for  metals  over  a 
strain  rate  range  of  500  -  2000/sec.  This  apparatus  is  also  useful  to  determine  the  compressive 
strength  of  unconfined  ceramics  under  uniaxial  stress  condition. 

2.1.1  Metals 

The  SHB  test  data  provide  variation  of  strength  with  respect  to  strain  rate  under  uniaxial 
stress.  Appendix  A  provides  the  SHB  data  for  the  various  metals  tested  under  this  program. 
Using  these  data,  the  Bodner-Partom  viscoplastic  constitutive  (strength)  model  [13]  constants 
were  determined.  Detailed  discussion  of  the  modeling  is  given  in  Section  3. 

The  IMACON  high  speed  camera  was  employed  to  obtain  data  for  ductile  necking 
evolution  in  copper,  Maraging  250  steel,  and  pure  tantalum.  Rajendran  and  Bless  [15]  described 
the  use  of  such  data  to  correct  the  engineering  stress-strain  relationship  beyond  the  onset  of 
necking.  The  measured  radius  of  curvature  of  the  neck  profile  and  the  minimum  radius  from  the 
high  speed  photographs  are  used  to  calculate  the  Bridgman  correction  factors.  Using  the 
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Bridgman  relationship,  the  SHB  stress-strain  curves  beyond  the  onset  of  necking  can  be 
constructed. 

2.1.2  Ceramics  and  Concrete 

The  compressive  strength  of  brittle  ceramics  and  concrete  materials  is  an  important 
mechanical  property  for  the  design  of  impact  resistant  structures.  Compressive  strength  data  can 
be  obtained  from  the  SHB  tests.  The  data  from  the  SHB  tests  are  also  useful  for  evaluating  the 
effects  of  strain  rate  on  the  strength  of  unconfined  brittle  materials.  We  performed  SHB  test' 
on  both  silicon  nitride  and  concrete.  The  results  for  concrete  which  were  reported  previously  by 
Antoun  [16]  and  John  et  al.  [17]  are  summarized  in  this  section.  The  compressive  strength 
variation  with  respect  to  strain  rate  for  silicon  nitride  and  concrete  was  determined  and  the  results 
are  presented  in  Figures  2.1  and  2.2,  respectively.  Both  the  ceramic  and  concrete  exhibited 
significant  rate  dependency.  For  concrete,  data  from  other  sources  are  also  included  for 
comparison.  The  data  shown  in  Figures  2.1  and  2.2  are  important  in  the  constitutive  model 
development. 

The  tensile  strength  of  brittle  materials  can  be  determined  from  ‘splitting  tension’  tests 
using  the  test  specimen  shown  in  Figure  2.3.  Quasi-static  and  dynamic  tests  were  conducted 
using  a  servo-hydraulic  testing  machine  and  a  split  Hopkinson  bar  (SHB),  respectively.  A 
schematic  of  the  SHB  setup  is  shown  in  Figure  2.4.  Using  the  equation  in  Figure  2.3,  the 
splitting  tensile  strength  is  calculated  at  the  peak  load  as  measured  by  the  strain  gauge  on  the 
SHB’s  transmitter  bar.  The  strain  rate  in  the  specimen  is  computed  using  the  peak  stress. 
Young’s  modulus,  and  time  to  peak,  assuming  elastic  behavior  [18].  The  compressive  loading 
results  in  a  uniform  lateral  (cx)  tensile  stress  in  the  middle  of  the  specimen  along  the  loading 
axis.  In  tension-weak  brittle  materials  such  as  ceramic  and  concrete,  tensile  failure  can  be 
expected  to  occur  with  the  crack  initiating  near  the  center  of  the  disk.  The  data  obtained  for 
silicon  nitride  are  plotted  in  Figure  2.5.  As  one  can  see  in  this  figure,  the  tensile  strength  is 
highly  strain  rate  sensitive. 

The  concrete  used  in  this  study  consisted  of  Type  I  ordinary  portland  cement,  river  gravel 
(maximum  size  =  0.16  in.)  and  water  in  the  ratio  1:2.5:0.5  by  weight  The  3  in.  diameter  x  6  in. 
cylinders  were  cast  in  cardboard  molds  and  cured  for  14  days  in  water  followed  by  14  days  in 
laboratory  air.  The  splitting  tension  specimens  were  prepared  from  some  of  these  cylinders. 
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Figure  2.1.  Silicon  Nitride  Compressive  Strength  Variation  with  Respect  to  Strain  Rate. 
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Figure  2.3.  Schematic  of  Splitting  Tension  Geometry. 
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Figure  2.5. 


Disk  specimens  with  diameters  of  0.5,  1.0  and  2.0  in.  and  thicknesses  of  0.25  and  0.50  in.  were 
prepared.  The  other  3x6  in.  cylinders  were  used  for  the  quasi-static  compression  tests.  The 
average  compressive  strength.  Young’s  modulus  (E),  and  Poisson’s  ratio  are  equal  to  8357  psi, 
4.25xl06  psi  and  0.21,  respectively. 

For  the  elastic  verification  tests,  a  2.0  in.  diameter  x  0.5  in.  thick  concrete  specimen  was 
instrumented  with  strain  gauges  at  the  center  on  both  sides.  This  specimen  was  first  loaded 
within  the  elastic  limit  at  the  quasi-static  rate  to  determine  the  Young’s  modulus  and  Poisson’s 
ratio.  The  same  specimen  was  then  subjected  to  low  velocity  impact  in  the  SHB.  The  tensile 
strain  measured  during  the  impact  event  is  shown  as  a  dotted  line  in  Figure  2.6.  Assuming  that 
the  deformation  was  elastic,  the  tensile  strain  at  the  center  of  the  disk  was  predicted  from  the 
output  of  the  gauge  on  the  SHB.  The  predicted  strain  was  then  time-shifted  and  plotted  as  a 
solid  line  in  Figure  2.6.  The  gauge  measurements  on  the  SHB  correlate  well  with  the  average 
specimen  response.  The  oscillations  in  the  strain  response  are  due  to  the  wave  reflections  within 
the  specimen.  The  average  time  period  of  the  oscillation  is  30  ps  which  is  equal  to  the  time 
taken  by  the  wave  to  traverse  twice  the  diameter  of  the  specimen,  assuming  that  the  Young’s 
modulus  of  concrete  is  rate-independent. 

In  the  disk  shown  in  Figure  2.3,  there  are  two  possible  locations  of  crack  initiation, 
namely  (1)  the  point  of  load  application  due  to  high  compressive  stresses  and  (2)  at  or  near  the 
center  of  the  disk  due  to  tensile  stresses.  To  determine  the  location  of  crack  initiation,  high 
speed  photography  tests  were  conducted.  The  concrete  specimens  were  painted  black  to 
emphasize  the  cracks.  During  the  impact,  the  photographs  were  taken  at  1.0  ps  intervals.  Two 
frames  highlighting  the  crack  initiation  event  and  the  corresponding  stress  versus  time  plot  are 
shown  in  Figures  2.7  and  2.8,  respectively.  The  first  frame,  at  a  time  equal  to  32  ps,  shows  that 
crack  initiation  occurred  in  the  constant  tensile  stress  region.  Hence  this  is  a  valid  splitting 
tension  test.  The  average  crack  velocity  calculated  from  these  tests  was  880  ft/s  (=0.12  x 
Rayleigh  wave  speed). 

The  results  of  the  quasi-static  and  dynamic  tests  are  plotted  in  Figure  2.9  [16,17].  The 
splitting  tensile  strength  of  concrete  is  shown  to  be  strain-rate  sensitive  as  was  also  observed  by 
Ross  et  al.  [18].  The  maximum  ratio  of  the  dynamic  to  quasi-static  strength  is  4.5  at  a  strain  rate 
of  73  per  second.  As  shown  in  Figure  2.9,  the  strain  rate  sensitivity  is  independent  of  thickness. 
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Figure  2.7. 
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Photographs  Highlighting  Locations  of  Crack  Initiation  in  Concrete. 
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ure  2.8.  Stress  History  Corresponding  to  the  Test  in  Figure  2.7  for  Concrete. 
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Figure  2.9.  Effect  of  Strain  Rate  and  Specimen  Size  on  the  Tensile  Strength  of  Concrete. 


Tensile  Strength  (MPa) 


The  quasi-static  strength  decreased  with  increasing  diameters  as  also  shown  by  other  investigators 
[18,19].  The  tensile  strength  from  the  2.0  in.  diameter  specimen  was  equal  to  about  one-tenth 
the  compressive  strength. 

The  dynamic  tensile  strength  data  was  normalized  with  respect  to  the  quasi-static  strength 
for  the  same  size  specimen  and  replotted  with  the  results  of  other  investigators  [18-21]  in  Figure 
2.10.  The  data  from  our  study  confirm  the  steep  increase  in  dynamic  tensile  strength  in  the 
strain  rate  regime  from  1.0  to  100  per  second. 

In  summary,  dynamic  splitting  tension  tests  were  conducted  using  the  SHB  apparatus  to 
investigate  the  influence  of  specimen  size  on  the  dynamic  tensile  st/ength  of  concrete.  The 
validity  of  the  SHB  data  was  confirmed  through  high  speed  photographs  of  the  concrete  specimen 
under  dynamic  deformation.  The  thickness  of  the  concrete  specimens  had  negligible  effect  on 
the  splitting  tensile  strength.  The  strength  decreased  with  increasing  specimen  diameter  under 
quasi-static  and  dynamic  loading.  But  the  ratio  of  dynamic  to  quasi-static  strength  was  observed 
to  be  size-independent.  For  strain  rates  greater  than  1.0  per  second,  the  splitting  tensile  strength 
exhibited  significant  strain-rate  dependency.  At  a  strain  rate  of  73  per  second,  the  ratio  of 
dynamic  to  quasi-static  strength  was  about  4.5. 

2.2  Plate  Impact  Experiments 

The  plate  impact  experiment  is  the  most  commonly  used  configuration  for  studying 
dynamic  tensile  (spall)  failure  in  materials  at  very  high  strain  rates.  In  the  plate  impact  test,  a 
flat  flyer  plate  is  made  to  impact  against  a  target  plate  at  a  high  velocity.  The  diagnostic 
measurements  usually  involve  a  VISAR  to  measure  the  particle  velocity  history  or  piezo-resistive 
stress  gauges  to  measure  the  stress  history  at  the  target  and  back  plate  interface.  This 
configuration  is  schematically  shown  in  Appendix  A.  The  flyer  and  target  may  be  of  the  same 
or  different  material.  Compressive  stresses  are  produced  and  transmitted  immediately  from  the 
plane  of  impact  to  the  adjacent  stress  free  areas  of  the  material  in  the  form  of  a  stress  pulse.  A 
typical  free  surface  velocity  history  is  shown  in  Figure  2.1 1.  The  portion  of  the  velocity  history 
indicated  by  the  letter  E  represents  the  arrival  of  an  elastic  shock  wave  from  the  impact  plane. 
The  material  particles  are  compressed  elastically  ( E-E 1 )  until  the  relatively  slow  moving  plastic 
shock  wave  arrives  at  the  free  surface.  Typically,  both  the  shape  and  time  duration  of  the  portion 
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DYNAMIC  STRENGTH  /  STATIC  STRENGTH 


STRAIN  RATE  (sec  ’) 


Figure  2.10.  Strain  Rate  Sensitivity  of  Tensile  Strength  Normalized  With  Respect  to  Static 
Strength  for  Concrete 
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Figure  2.11.  Typical  Free  Surface  Velocity  of  the  Target  Measured  Using  VISAR 


E-E /  depend  on  the  strain  rate  sensitivity  of  the  material  and  thickness  of  the  target  plate.  For 
engineering  calculations,  the  average  velocity  level  between  the  points  E  and  E1  is  usually 
defined  as  the  velocity  at  the  HEL,  Vp^.  Using  this  experimental  value,  the  stress  at  the  HEL 
(°HEl)  *s  calculated  from  the  relationship  =  '/Spc V^j ,  where  p  is  the  material  density,  c 

is  the  elastic  sound  wave  speed,  and  is  the  velocity  at  the  HEL  measured  at  the  target’s 
free  surface.  This  expression  is  valid  only  under  symmetric  impact  conditions,  that  is,  when  the 
target  and  flyer  are  of  the  same  material.  The  high  strain  rate  yield  strength,  Y0,  can  then  be 

K  2 

calculated  as  Y0  =  Ohel  /  [  +  ^.  ],  where  G  is  the  shear  modulus  and  K  is  the  bulk 

modulus. 

When  the  plastic  shock  wave  arrives  at  the  target’s  free  surface,  this  wave  takes  the 
particle  velocity  to  the  peak  level  defined  between  points  P  and  R.  Later,  the  elastic  release 
wave  from  the  back  surface  of  the  flyer  plate  arrives  at  point  R,  and  the  velocity  drops  to  R 1 . 
The  plastic  release  wave  from  the  flyer  reaches  the  target’s  free  surface  and  drops  the  velocity 
further  to  F.  In  general,  in  modeling,  the  level  E-E '  and  the  wave  structure  around  the  HEL 
are  sensitive  to  the  constitutive  relations  (strength  models);  however,  the  slopes  and  peak  of  the 
VISAR  signal  depend  primarily  on  the  pressure-volume  relationship  (equation  of  state). 

Depending  on  the  severity  of  the  impact,  the  bulk  material  (target)  can  spall  and  a  stress- 
free  fracture  surface  (spall  plane)  is  created  in  the  target.  When  the  spall  plane  separates  inside 
the  target,  the  stresses  relax  to  zero,  and  hence,  a  compressive  shock  wave  is  generated.  This 
shock  wave  arrives  at  point  S,  as  shown  in  the  typical  VISAR  signal.  The  material  particle 
velocity  history  follows  the  dashed  line  (recompression).  The  signal  beyond  point  S  is  often 
denoted  as  the  "spall  signal."  This  spall  signal  is  most  often  used  in  the  estimation  or  calibration 
of  failure  model  parameters. 

The  experimental  data  from  the  plate  impact  experiments  for  several  metals  and  ceramics 
are  given  in  Appendix  A. 
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2.2.1  Metals 


We  performed  plate  impact  experiments  on  several  steels,  aluminum  alloys,  copper,  pure 
tantalum,  and  tungsten.  The  details  of  the  tests  are  given  in  Table  2.1.  The  stress  gauge 
measured  stress  histories  at  the  interface  of  the  target  and  the  back  plate  are  given  in  Appendix 
A.  The  measured  stress  history  is  often  used  in  the  spall  model  calibration.  The  Hugoniot  elastic 
limits  (HEL)  for  these  various  metals  were  determined  and  tabulated  in  Table  2.2. 


TABLE  2.1 

PLATE  IMPACT  EXPERIMENTS  ON  METALS 


Shot 

Number 

Flyer  Plate 

Target  Plate 

Impact 

Velocity 

(m/sec) 

Material 

Thickness 

(mm) 

Material 

Thickness 

(mm) 

7-538 

Copper 

2.0 

Copper  (98%) 

9.0 

185 

7-1070 

Copper 

2.0 

Tungsten  (98%) 

4.0 

188 

7-1071 

46100  Steel 

2.65 

46100  Steel 

6.0 

310 

7-1244 

Copper 

6.4 

Tungsten  (90%) 

5.6 

292 

7-1250 

2519  A1 

4.0 

2519  A1 

1 

324 

7-1267 

Copper 

4.0 

1215  Steel 

M 

634 

7-1268 

1020  Steel 

3.9 

1020  Steel 

7.8 

572 

7-1288 

Copper 

2.0 

Copper 

4.1 

512 

7-1298 

Copper 

3.0 

Armco  Iron 

6.0 

470 

7-1299 

Copper 

3.0 

HY100  Steel 

6.0 

489 

7-1300 

Copper 

3.0 

Cl 008  Steel 

6.0 

462 

7-1302 

Copper 

2.0 

Tantalum 

6.0 

735 

7-1305 

Copper 

3.0 

4340  Steel 

665 

7-1454 

MAR-250  Steel 

3.0 

MAR-250  Steel 

i 

575 

7-1455 

MAR-200  Steel 

3.0 

MAR-200  Steel 

6.0 

508 

7-1523 

AF1410  Steel 

mm 

AF1410  Steel 

9.0 

678 
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TABLE  2.2 
HEL  FOR  METALS 


: 


MATERIAL 

HEL  (kbar) 

Y  (kbar) 

Tantalum 

27.0 

13.5 

Tungsten 

82.0 

51.0 

46100  Steel 

38.0 

19.0 

1215  Steel 

14.4 

7.2 

4340  Steel 

20.0 

10.0 

MAR-200  Steel 

26.0 

13.0 

MAR-250  Steel 

21.0 

10.5 

1020  Steel 

7.4 

3.7 

C1008  Steel 

11.0 

5.5 

AF1410  Steel 

36.0 

18.0 

HY100  Steel 

16.0 

8.0 

2.2.2  Ceramics 

We  performed  a  limited  number  of  plate  impact  tests  on  AD85,  AD90,  AD95,  pure 
alumina  (99.8%),  Silicon  Nitride,  Boron  Carbide,  and  Titanium  Diboride.  The  details  of  the 
experiments  are  given  in  Table  2.3.  The  HEL  values  for  various  ceramics  are  summarized  in 
Table  2.4.  The  corresponding  measured  stress  histories  are  given  in  Appendix  A. 

2.3  Plate-on-Bar  Impact  Data 

In  this  configuration,  a  thick  flyer  plate  impacts  onto  a  long  slender  bar.  An  aspect  ratio 
(length  to  diameter)  of  twelve  was  used  for  the  bar.  In  Figure  2.12,  a  schematic  of  the  plate-on- 
bar  test  is  shown.  For  diagnostic  measurement,  a  piezo-resistive  stress  gauge  is  imbedded  into 
the  bar  at  a  distance  several  diameters  away  from  the  impact  end.  This  technique  was  described 
in  detail  by  Rosenberg  and  Bless  [22]. 


TABLE  2.3 

PLATE  IMPACT  EXPERIMENTS  ON  CERAMICS 


Shot 

Number 

Flyer  Plate 

Target  Plate 

Impact 

Velocity 

(m/sec) 

Material 

Thickness 

(mm) 

Material 

Thickness 

(mm) 

7-1108 

Copper 

2.6 

AD-90 

I 

i 

I 

181 

7-1109 

Copper 

2.6 

AD-99.5 

8.4 

551 

7-1136 

Copper 

2.6 

TiB2 

12.7 

719 

7-1154 

Copper 

2.0 

B4C 

10.0 

1121 

7-1215 

Copper 

3.0 

Si3N4 

11.4 

1009 

7-1218 

Copper 

3.0 

Si3N4 

11.4 

1110 

7-1295 

AD-85 

5.0 

AD-85 

8.0 

237 

7-1332 

Copper 

3.0 

Si3N4 

11.4 

770 

TABLE  2.4 
HEL  FOR  CERAMICS 


MATERIAL 

HEL  (kbar) 

Y  (kbar) 

Alumina  (AD-85) 

55 

39 

Alumina  (AD-90) 

70 

52 

Alumina  (AD-99.5) 

83 

59 

Alumina  (Hot  Pressed,  99.8%  Pure) 

123 

82 

Silicon  Nitride 

120 

78 

Silicon  Carbide 

140 

107 

Aluminum  Nitride 

94 

64 

Titanium  Diboride 

75 

66 

Boron  Carbide 

180 

146 

Sodalime  Glass 

64 

39 
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2.3.1  Metal  Bars 


While  the  stress  state  near  the  impact  end  is  multiaxial,  it  is  uniaxial  at  the  gauge 
location.  At  high  velocities,  metals  plastically  flow  near  the  impact  end.  However,  the  stress 
wave  that  passes  through  the  gauge  is  essentially  an  elastic  wave  and  its  amplitude  is  a  measure 
of  yield  strength  at  medium  strain  rates  (about  500/sec).  This  measurement  is  useful,  especially 
when  the  SHB  test  does  not  provide  accurate  measurements  of  the  initial  yield  strength  due  to 
spurious  wave  reflections  at  small  strains  (<  0.03). 

Grove  and  Rajendran  [23]  performed  a  numerical  simulation  of  a  15mm  thick  1020  steel 
plate  impacting  a  152mm  long,  12.7mm  diameter  Cl 008  steel  bar.  The  calculated  initial  yield 
stress  (maximum  amplitude  of  the  elastic  wave)  matched  well  with  the  gauge  data.  Such 
numerical  simulations  are  useful  in  the  evaluation  of  constitutive  models. 

2.3.2  Ceramic  Bars 

The  plate-on-bar  experimental  technique  was  initially  employed  to  measure  stresses  in 
metal  bars.  However,  under  the  Air  Force  contract,  we  extended  the  technique  to  characterize 
deformation  and  fracture  in  ceramic  bars.  Several  articles  were  prepared  and  published  by  Brar 
et  al.  [24],  Rosenberg  et  al.  [25],  Brar  and  Bless  [26],  Grove  and  Rajendran  [27],  and  others 
[28,29]. 

Brar  et  al.  [24]  conducted  plate-on-bar  experiments  on  Coor’s  AD-998  (almost  pure 
aluminum  oxide)  and  AD94.  Ten  experiments,  at  impact  velocities  in  the  range  96-560  m/s,  were 
performed  using  a  50  mm  gas/powder  gun.  Steel  impactors,  15  mm  thick,  were  used  in  all  the 
experiments.  In  these  experiments,  when  the  impact  stress  was  below  the  compressive  yield 
strength,  a  constant  amplitude  stress  wave  was  produced  in  the  bar.  When  the  impact  stress 
exceeded  the  yield  strength,  the  stress  level  in  the  bar  decayed  with  distance  and  time.  This 
stress  decay  is  consistent  with  lateral  release  involving  dilatancy. 

These  types  of  impact  tests  can  be  used  to  measure  dynamic  yield  strength.  The 
measured  dynamic  strength  in  AD998  and  AD94  ceramics  did  not  differ  significantly  from  static 
strengths.  Typical  oscilloscope  records  from  experiments  at  impact  velocities  of  102,  125,  and 
551  m/s  respectively  on  AD998  (99.8  %  pure  alumina)  bars  are  shown  in  Figure  2.13.  The 
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Figure  2.13.  Oscilloscope  Records  from  Experiments  on  AD-998  Rods.  Scales  for  All  Records 
are  0.2  V/div.  and  0.5  ps/div. 
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corresponding  peak  stresses  were  18,  21.5,  and  92  kbars.  Detailed  interpretations  of  these  plate- 
on-bar  experiments  can  be  found  in  Reference  [24]. 

The  serial  high  speed  photographs  of  the  fracturing  AD998  bar  are  shown  in  Figure  2.14. 
The  frames  are  read  from  top  to  bottom  and  left  to  right.  The  first  frame  (top  left)  is  at  time 
zero  and  successive  frames  are  taken  every  10  microseconds.  In  the  experiment,  the  bar  was 
painted  black  to  observe  the  fracture  patterns.  The  bar  laterally  expands  due  to  pulverization  at 
the  impact  end.  Later,  splitting  cracks  emanate  from  the  impact  end  and  propagate  towards  the 
other  end  of  the  bar.  The  fracture  of  the  bar  is  indeed  extremely  complex  due  to  fracturing  under 
a  multiaxial  stress  state.  At  present,  no  definitive  understanding  or  modeling  of  this  failure 
process  is  available. 

2.4  Bar-on-Bar  Impact  Test 

Failure  of  ceramics  can  be  studied  with  instrumented  bar  impacts  and  high  speed 
photography.  Brar  and  Bless  [26]  summarized  their  results  in  detail.  In  the  bar-on  bar  impact 
configuration  (see  Figure  2.15),  a  short  ceramic  bar  impacts  onto  a  long  ceramic  bar  having  an 
aspect  ratio  of  4~5.  Both  bars  are  made  of  the  same  material.  The  strain  rates  are  in  the  range 
of  103  ~  lO'Vsec.  Manganin  gauges  were  embedded  at  several  diameters  away  from  the  impact 
end.  The  rod  projectiles  were  launched  using  a  lexan  sabot  in  a  50  mm  gas/powder  gun.  The 
bar  targets  were  aligned  for  a  planar  impact  using  a  special  fixture.  An  Imacon  790  high  speed 
camera  was  used  to  photograph  the  fracture  of  the  rods.  Alumina  bars  were  painted  black  so  that 
the  cracks  and  faults  could  be  distinguished.  The  bar  fracturing  was  photographed  every  10 
microseconds.  The  recording  time  duration  is  usually  about  140  microseconds.  The  photographs 
of  the  fracturing  ceramics  indicate  three  types  of  fracture;  pulverization  of  the  bar  at  the  impact 
plane,  axial  splitting  (or  faulting),  and  a  propagating  destruction  (fracture)  wave.  The  material 
behind  this  wave  is  pulverized  and  the  measurements  from  the  pictures  indicate  that  the  fracture 
wave  propagates  faster  than  the  interface  velocity. 

In  alumina  AD-998,  a  splitting  type  of  crack  pattern  can  be  seen  from  Figure  2.16. 
However,  in  soda  lime  glass  bars  the  splitting  type  fracture  was  not  observed  in  the  high  speed 
photographs  shown  in  Figure  2.17.  A  compressive  fracture  front  propagates  from  the  impact  end. 
This  fracture  pattern  is  more  of  a  pulverization  type  than  a  few  cracks  running  into  the  bar  as 
with  splitting  cracks.  Upon  impact,  a  compressive  shock  wave  travels  into  the  bar  and  reflects 
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Figure  2. IS.  A  Schematic  of  Bar- 
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Figure  2.17.  High  Speed  Photograph  of  Soda  Lime  Glass  Fracturing  Due  to  Soda  Lime  Glass 
Bar-on-Bar  Impact. 
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as  a  tensile  wave.  This  tensile  wave,  along  with  release  waves  from  the  lateral  surface,  initiates 
fracture  at  the  free  end  of  the  bar.  The  tension-initiated  fracture  front  propagates  towards  the 
fracture  front  arriving  from  the  impact  end  of  the  bar. 

This  tensile  wave,  which  propagates  faster  than  the  fracture  front,  encounters  the 
compressive  fracture  front  and  reflects  back  as  a  compressive  wave  of  lesser  amplitude.  This 
trapped  wave  decays  and,  as  one  can  clearly  see  in  Figure  2.17,  leaves  a  small  length  of  the  bar 
unbroken.  These  types  of  diagnostic  photographs  should  be  useful  to  model  builders  in  validating 
ceramic  failure  models. 
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Section  3 

Ductile  Failure  Model  for  Metals 


Three  of  the  main  features  in  a  dynamic  ductile  failure  model  are  the  following:  (1)  the 
initial  intact  (void-free)  material  requires  a  constitutive  description  which  will  include  strain  rate 
and  temperature  effects;  (2)  a  mathematical  description  of  the  void  nucleation  and  growth 
process;  and  (3)  pressure  dependent  plastic  flow  equations  for  the  porous  (void-containing) 
aggregate.  There  could  be  several  other  features  in  the  failure  modeling.  For  instance,  depending 
on  the  level  of  voids,  the  stiffness  of  the  aggregate  may  significantly  be  reduced;  therefore, 
models  to  degrade  the  stiffness  will  be  required  for  the  realistic  description  of  the  material 
behavior.  The  final  process  of  voids  coalescence  may  also  be  important.  Another  aspect  of  the 
failure  model  is  the  equation  of  state  for  the  porous  aggregate. 

In  this  section,  a  brief  description  of  a  ductile  failure  model  developed  under  this 
program  is  presented.  Previously,  a  detailed  description  of  this  RDG  model  was  presented  by 
Rajendran  [8],  and  Rajendran,  Dietenberger,  and  Grove  [1].  This  three-dimensional,  continuum- 
mechanics  based,  dynamic  failure  model  (RDG  Model)  is  capable  of  describing  spallation  under 
two-dimensional  stress  (long  rod  impact  on  a  thick  plate)  and  one-dimensional  strain  (plane  plate 
impact).  The  RDG  model  considered  a  viscoplastic  constitutive  description  for  the  matrix  and 
the  porous  aggregate  materials.  The  stress  and  strain  based  void  nucleation  process  was  modeled 
through  a  Gaussian  distribution.  The  RDG  model  proposed  a  new  pressure  dependent  yield 
function  for  describing  plastic  flow  in  the  porous  aggregate.  There  are  four  phases  in  the  model. 
In  the  first  phase,  the  intact  material  is  described  by  the  Bodner-Partom  viscoplastic  model  [13]. 
The  void  nucleation  is  introduced  in  the  second  phase.  The  void-contained  aggregate  is  described 
in  the  third  phase  using  an  associated  plastic  flow  rule  derived  from  a  pressure  dependent  yield 
function.  The  last  phase  of  modeling  is  the  coalescence  of  voids  leading  to  complete  failure. 
In  the  RDG  model,  separate  modeling  of  the  coalescence  process  is  not  needed.  The  void  growth 
law  is  such  that  the  growth  rate  is  rapidly  increased  as  the  damage  approaches  its  critical  value. 
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3.1  Constitutive  Equations 


A  pressure  dependent  yield-criterion-based  approach  has  been  considered  in  the 
constitutive  model  formulation.  For  randomly  distributed  voids  or  microcracks  contained  in  the 
aggregate,  the  yield  behavior  will  be  influenced  not  only  by  the  second  invariant  of  the  deviatoric 
stress  (J2)  but  also  by  the  pressure  or  mean  stress  (Ij).  The  following  form  of  the  yield  function 
has  been  considered: 


<I>=  (2+p2)  J2 


(1) 


where. 


5  (p)  = 


g(p)  -g(pcr) 
g(D  -  g(pcr) 


(2) 


and 


g(p)’(K_w[(i'p)r  <3) 

where  Ym  is  the  effective  yield  stress  in  the  matrix  material,  K  and  N  are  model  constants,  and 
p  is  the  relative  density.  A  negative  value  of  N  makes  Equation  (3)  a  hyperbolic  power  function. 
Numerical  simulations  of  a  plate  impact  test  configuration  were  performed  to  evaluate  the  effects 
of  the  8  function  on  the  spall  signal.  In  these  simulations,  the  initial  slope  of  8(p)  at  p=l  greatly 
influenced  the  slope  of  the  spall  signal.  Consequently,  the  parameter  P  (=5/  (1))  was  introduced. 
With  p,  N,  and  pcr  as  model  constants,  the  corresponding  value  of  K  can  be  solved  by  a  simple 
iterative  scheme.  An  idealistic  value  of  zero  can  be  assumed  for  pcr. 

The  viscoplastic  strain  rates,  in  the  aggregate  can  be  calculated  using  the  flow  rule 
derived  from  the  yield  function  given  by  Equation  (1).  An  expression  for  Gjj  can  be  obtained 
from  Hooke’s  law  (elastic  stress-strain  relationship)  by  defining  the  elastic  strain  rate  as  the 
difference  between  total  and  plastic  strain  rates. 


G  ij  ~  E  ik  kj  ~  ^  /c j ) 


(4) 


3-2 


The  evolution  laws  for  the  void  content  f  (void  volume  fraction,  f  =  1-p)  are  given  by 
the  nucleation  and  growth  equations.  Rajendran  et  al.  [1],  described  the  evolution  laws  in  detail. 
For  completeness,  the  salient  equations  are  given: 


where,  for  the  void  volume  fraction  nucleation  law  we  use 

=  (Ym  +  P)  +FtD% 


with 


(5) 


(6) 


(7) 


and 


f2 

s2yf2% 


(8) 


where  Ym  is  the  effective  stress  in  the  matrix,  P  is  the  pressure  in  the  aggregate,  DP  is  the 
equivalent  plastic  strain  in  the  matrix,  and  aN  and  eN  are  the  mean  equivalent  stress  and  strain, 
respectively,  around  which  the  nucleation  stress  and  strain  are  distributed  in  a  Gaussian  manner. 
The  terms  Sj  and  s2  are  the  standard  deviations  of  these  distributions.  These  two  parameters 
control  the  ranges  of  stress  or  strain  over  which  most  of  the  voids  can  be  nucleated.  The  terms 
fj  and  f2  define  the  maximum  allowable  void  volume  fractions  due  to  stress  and  strain  nucleation, 
respectively. 

The  growth  law  can  be  directly  related  to  the  dilatation  due  to  growth  of  voids  in  the 
aggregate.  By  definition,  the  void  growth  rate  is  given  by, 

fg  =  ptPU  <9> 
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where  repeated  index  means  summation,  6^  are  the  plastic  strain  rates  in  the  three  principal 
directions,  and  p  (=l-f)  is  the  relative  density. 

To  complete  the  formulation,  an  evolution  equation  for  the  matrix  effective  stress  is 
required.  For  the  intact  matrix  material,  either  the  Bodner-Partom  [13]  (BP)  or  Johnson-Cook 
[14]  (JC)  model  can  be  employed.  These  models  are  described  in  Appendix  B.  A  strain  rate 
dependent  relationship  between  Ym  and  the  equivalent  plastic  strain  rate  in  the  matrix  material 
is  provided  by  these  models. 

3.2  Degradation  of  Shear  and  Bulk  Moduli 

The  shear  and  bulk  moduli  (K  and  G,  respectively)  are  degraded  using  Mackenzie’s 
formula  [30].  Mackenzie  derived  an  approximate  analytical  expression  for  the  reduction  of 
elastic  stiffness  based  on  the  elastic-plastic  flow  around  a  spherical  void  in  an  infinite, 
incompressible  matrix  material.  The  RDG  model  [1]  employed  the  corresponding  expression  for 
the  degraded  moduli  K  and  G  as. 


K  _ 
~K 


(1  -f) 


3K^> 


(1 


(10) 


and. 


(1-f) 


J-  _  6K  +  12G  ^1 
\  9k+8g  ) 


(11) 


where  K  and  G  are  the  bulk  and  shear  moduli  of  the  intact  material.  Appendices  C  and  D  outline 
the  solution  scheme  to  solve  the  governing  equations. 

3.3  Equation  of  State  for  the  Aggregate 

The  equation  of  state  for  the  solid  (intact)  material  has  been  routinely  determined  from 
the  experimentally  obtained  Hugoniot  data  for  steel,  aluminum,  copper,  and  other  metals. 
However,  once  the  initially  intact  solid  develops  voids  upon  tensile  loading,  the  equation  of  state 
(EOS)  of  the  original  solid  material  is  no  longer  valid  during  the  subsequent  compressive  loading. 
For  this  purpose,  the  EOS  parameters  have  to  be  modified. 
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The  simple  approach  used  in  the  RDG  model  is  an  extension  of  Mackenzie’s  procedure. 
The  EOS  for  the  intact  solid  is  described  by  the  Mie-Gruneisen  equation  as 

ps=  (p14+p2^2+p3^i3)  a-r£)  +r(j-i0)  (i2) 

where  Pj,  P2,  P3,  and  T  are  EOS  constants,  p  is  the  volumetric  compressible  strain  of  the  intact 
solid,  (  - 1 )  ,  and  I  is  the  specific  internal  energy  of  the  solid  material.  To  account  for  the 

voids,  Mackenzie’s  formulation  is  implemented  as  follows.  A  volumetric  elastic  strain  of  the 
aggregate  is  first  defined  as, 

(13) 

to  replace  p  in  the  right  side  of  Equation  (12).  Then  the  right  side  expression  of  Equation  (12) 

is  multiplied  by  the  Mackenzie  correction  term,  —  ,  to  obtain  the  aggregate  pressure,  P.  Thus 

K 

a  voided  material  undergoing  compression  will  load  along  a  degraded  bulk  modulus  at  the 
aggregate  compressible  strain.  At  high  enough  stress  levels,  the  void  will  collapse  (f=0),  and  thus 
the  modified  EOS  will  reduce  to  Equation  (12). 

The  various  governing  equations  can  be  rearranged  and  numerically  integrated.  The  RDG 
model  along  with  the  choice  of  BP  or  JC  model  was  implemented  into  the  EPIC-2  finite  element 
code  [31].  Special  purpose  subroutines  have  been  developed  successfully.  The  corresponding 
numerical  scheme  is  given  in  Appendices  C  and  D. 

The  strength  model  (BP  or  JC)  constants  are  determined  using  the  stress  -  strain  data 
from  split  Hopkinson  bar  (SHB)  tests.  These  data  may  comprise  data  from  quasi-static  tensile, 
SHB  tensile,  compressive,  and  torsional  tests.  Rajendran  et  al.  [32]  described  a  combined 
experimental  and  numerical  scheme  for  evaluating  BP  model  constants.  They  also  included  the 
HEL  (Hugoniot  Elastic  Limit)  data  in  the  scheme.  Johnson  and  Holmquist  [33]  provided  a 
methodology  to  obtain  JC  model  constants.  The  BP  and  JC  model  constants  for  various  metals 
are  given  in  Appendix  B. 
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Determining  failure  model  constants  requires  plate  impact  test  data  in  terms  of  either 
velocity  history  or  stress  history  at  the  back  of  the  target.  Rajendran  et  al.  [1]  determined  the 
failure  model  constants  using  the  velocity  history  obtained  from  a  plate  impact  test  on  OFHC 
copper. 

3.4  Model  Parameters 

The  BP  model  was  employed  for  the  intact  material  description.  The  corresponding  model 
parameters  are  given  in  Appendix  B.  For  establishing  a  standard  procedure  to  calibrate  the 
failure  model  constants,  two  plate  impact  experiments  (7-538  and  7-1299)  were  considered.  (See 
Table  2.1.)  Experiment  7-538  was  first  employed  in  the  failure  model  parameter  evaluation. 
In  this  experiment,  a  2  mm  copper  flyer  was  impacted  against  a  9  mm  OFHC  copper  target  at 
an  impact  velocity  of  185  m/s,  and  the  target’s  free  surface  velocity  history  was  obtained  using 
a  VISAR  (velocity  interferometer).  The  second  experiment  (7-1299)  consisted  of  a  3  mm  copper 
flyer  impacting  a  6  mm  HY100  steel  target  at  a  velocity  of  489  m/s;  in  this  experiment,  the 
diagnostic  was  a  manganin  gauge  measured  stress  history.  These  two  experiments  were 
simulated  using  the  EPIC-2  code.  The  special  purpose  subroutines  describing  the  RDG  model 
were  used  to  model  the  spall  process  in  both  OFHC  copper  and  HY100  steel.  Using  the 
suggested  procedure,  the  model  constants  for  these  two  as  well  as  several  other  metals  were 
determined.  The  corresponding  constants  are  given  in  Table  3.1. 

The  comparison  between  the  model  simulation  and  the  experimental  velocity  history  data 
for  OFHC  copper  is  shown  in  Figure  3.1.  Since  the  HEL  of  OFHC  copper  is  extremely  low 
compared  to  the  shock  stress,  the  BP  model  predicted  HEL  (knee  at  A  in  Figure  3.1)  did  not 
match  well.  However,  for  copper,  the  HEL  value  is  not  considered  critical  because  plastic 
yielding  occurs  at  a  (uniaxial)  stress  below  1  kbar.  The  matching  beyond  point  S  is  controlled 
by  the  RDG  model  constants.  In  general,  the  overall  match  between  the  simulation  and  the  data 
is  very  good.  This  demonstrates  the  ability  to  calibrate  the  model  constants  using  plate  impact 
test  data  as  well  as  the  RDG  model  capability  to  accurately  reproduce  the  spall  signal. 


3-6 


TABLE  3.1 

RDG  MODEL  CONSTANTS 


Material 

Nucleation 

Yield  Function 

h 

aN 

(GPa) 

s(=aN/4) 

(GPa) 

P 

N 

OFHC  Copper 

0.01 

1.6 

0.4 

65 

-2.4 

HY100  Steel 

0.035 

3.6 

0.9 

30 

-1.0 

C1008  Steel 

0.02 

1.2 

0.3 

100 

-2.0 

MAR-200  Steel 

0.01 

4.8 

1.2 

40 

-5.0 

Armco  Iron 

0.05 

2.0 

0.5 

70 

-0.5 

Tantalum 

0.003 

4.8 

1.2 

10 

-3.0 

1020  Steel 

0.01 

2.4 

0.6 

1 

-0.5 

MAR-250  Steel 

0.01 

3.2 

0.8 

50 

■ 

o 

AF1410  Steel 

0.01 

5.6 

1.4 

20 

-4.0 

To  further  verify  the  model  constant  calibration  scheme  using  data  from  a  plate  impact 
test,  spall  in  HY100  steel  was  considered.  In  this  case,  a  manganin  gauge  recorded  stress  history 
provided  the  experimental  data.  Recall  that  the  failure  model  constants  for  OFHC  copper  were 
obtained  using  the  velocity  history.  The  experimental  data  were  obtained  only  for  a  time  duration 
of  4  microseconds.  (In  this  experiment,  release  waves  from  the  edges  of  the  flyer  plate  arrived 
at  the  gauge  location  at  around  4.2  microseconds.)  The  model  constants  were  obtained  by 
matching  the  experimental  stress  signal  until  this  time.  The  corresponding  match  is  shown  in 
Figure  3.2a.  The  excellent  comparison  between  model  and  experiment  before  the  spall  signal 
indicates  that  the  BP  model  constants  reproduced  the  HEL  data  well.  Also,  the  strain  rate 
sensitivity  of  HY100  steel  seemed  to  be  modeled  correctly  as  one  can  interpret  from  the  shape 
of  the  simulated  signal  around  the  stress  peak.  Rajendran  et  al.  [32]  modeled  the  SHB  test 
accurately  using  the  same  BP  model  constants  for  HY100  steel.  The  ability  to  model  both  the 
SHB  and  plate  impact  tests  indicates  the  generality  of  the  model  constants.  In  Figure  3.2b,  the 
computed  stress  history  at  the  spall  plane  is  shown.  The  compressive  stress  is  shown  as  positive, 
and  the  tensile  stress  is  shown  as  negative.  The  material  experiences  a  compressive  stress  of 
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Comparison  between  Model  and  Experimental  Stress  Histories  (HY100  Steel). 
Numerical  Simulation  of  the  Stress  History  at  the  Spall  Plane  (HY100  Steel). 


around  10  GPa.  When  the  material  goes  into  tension,  the  stress  peak  only  reaches  2.8  GPa.  This 
reduction  in  stress  under  tension  is  due  to  void  softening.  The  stress  relaxes  to  zero  as  the 
material  spalls.  Modeling  of  this  time  dependent  stress  relaxation  is  the  most  important  feature 
of  the  RDG  model. 

The  failure  model  constants  were  accurately  determined  as  can  be  seen  from  the  excellent 
comparison  between  the  model  simulation  and  the  experimental  data  in  Figure  3.2a.  Since  the 
BP  model  constants  for  the  intact  material  and  the  RDG  model  constants  for  the  void-containing 
aggregate  were  accurately  determined,  the  entire  stress  history  could  be  reproduced  accurately. 
Of  course,  this  statement  assumes  that  the  equations  of  state  for  the  flyer,  target,  and  the  PMMA 
are  accurately  modeled  in  the  computer  code  simulation.  The  ability  of  the  model  to  reproduce 
the  entire  spall  signal  clearly  demonstrates  that  the  nucleation  and  growth  models  follow  the 
physical  processes  realistically. 

RDG  model  constants  were  also  determined  for  1020,  C1008,  MAR-250,  AF1410,  and 
MAR-200  steels,  Armco  Iron,  and  pure  tantalum  by  matching  the  experimental  data  as  shown 
in  Figures  3.3  and  3.4.  The  ability  of  the  RDG  model  in  modeling  the  spall  failure  is  truly 
outstanding.  The  RDG  model  generated  spall  signal  for  tantalum  matched  extremely  well  with 
the  data  up  to  3.4  microseconds.  The  difference  between  the  model  and  experiment  could  be  due 
to  the  uncertainty  in  the  experimental  data.  In  this  test  (7-1302),  the  plate  impact  configuration 
was  such  that  the  spall  plane  location  was  very  close  to  the  manganin  gauge  (about  one 
millimeter). 

In  the  simulations,  behavior  of  the  intact  material  can  also  be  described  by  the  JC  model 
instead  of  the  BP  model.  Using  the  JC  and  RDG  models,  the  plate  impact  test  7-1299  on  HY100 
steel  was  resimulated.  The  JC  constants  for  the  HY100  steel  were  obtained  by  Johnson  and 
Holmquist  [33]  using  torsional  SHB  test  data.  Initially,  the  RDG  model  constants  used  were 
those  that  had  been  previously  determined  for  HY100  steel  in  conjunction  with  the  BP  model. 
Figure  3.5  compares  the  stress  histories  obtained  from  (1)  the  test  data,  (2)  the  simulation  using 
the  BP  model  with  the  RDG  model,  and  (3)  the  simulation  using  the  JC  model  with  the  RDG 
model.  In  the  simulation  using  the  JC  model,  the  spall  signals  did  not  match  exactly.  This  slight 
difference  could  be  attributed  to  the  material  characterizations  of  HY100  steel  using  different 
experimental  configurations  and  sources.  In  other  words,  the  descriptions  of  HY100  steel  by  the 
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Figure  3.5.  Comparison  of  Spall  Modeling  for  HY100  Steel  using  BP  and  JC  Models 
Spall  Process  is  Described  by  the  RDG  Model. 
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BP  model  and  the  JC  model  are  slightly  different;  because  of  this,  accurate  modeling  of  spallation 
using  the  JC  model  required  minor  adjustments  to  the  RDG  model  constants.  To  determine  the 
RDG  model  constants  that  correspond  to  the  use  of  the  JC  model  for  the  intact  material,  plate 
impact  simulations  were  made  and  the  constants  were  determined.  An  adjustment  to  the  fj  from 
0.035  to  a  value  of  0.01  led  to  a  better  fit  between  the  failure  model  and  the  test  data.  It  can  be 
seen  from  Figure  3.6  that  the  spall  signal  is  matched  extremely  well  by  both  BP  and  JC  models 
with  the  use  of  appropriate  RDG  model  constants. 

An  important  aspect  of  the  modeling  is  the  model’s  ability  to  reproduce  an  experiment 
which  was  not  used  for  determining  the  model  constants.  Also  the  model’s  ability  to  reproduce 
a  failure  process  in  some  other  stress-strain  configuration  is  important.  To  check  this  aspect,  a 
new  experiment  (7-1288)  on  OFHC  copper  was  considered.  In  this  experiment,  the  stress  history 
at  the  interface  of  the  target  and  the  back  plate  was  measured  using  a  manganin  gauge.  The 
target  thickness  and  the  impact  velocity  were  both  different  than  in  experiment  7-538.  To  check 
the  model’s  ability  to  predict  the  measured  stress  history,  the  new  experiment  7-1288  was 
simulated  using  the  same  model  constants  that  were  determined  using  experiment  7-538.  In 
Figure  3.7,  a  comparison  between  the  model  prediction  and  the  data  is  shown.  The  model 
reproduced  the  measured  spall  signal  extremely  well  as  can  be  seen  from  this  figure. 

3.5  Modeling  of  Double  Flyer  Plate  Impact 

In  application  problems,  dealing  with  projectile  penetration  into  a  target,  explosively 
compacted  metal  operations,  etc.,  the  loading  conditions  often  lead  to  multiple  shocking  in  the 
material.  Under  these  conditions,  both  void  growth  and  void  collapse  may  occur.  The  dynamic 
properties  of  the  shocked  material  are  modified  due  to  loading  history  effects.  To  simulate  such 
conditions  and  to  further  evaluate  the  RDG  model,  a  double  flyer  plate  impact  experiment  was 
simulated  with  the  EPIC-2  finite  element  code  [31].  Previously  calibrated  RDG  model 
parameters  were  used  in  the  simulation.  The  computed  stress  history  was  compared  with  the 
experimentally  measured  stress  history.  A  simple  void  collapse  criterion,  based  on  a  critical  void 
volume  fraction,  was  employed  to  obtain  a  better  fit  to  the  experimental  data. 

Yaziv  [34]  and  Yaziv  et  al.  [35]  introduced  double  flyer  plate  impact  techniques  for 
examining  the  dynamic  properties  of  shock  damaged  materials.  This  technique  differs  from  a 
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Figure  3.6.  Comparison  between  the  Spall  Signals  Generated  using  JC  and  BP  Models  for 
Different  Values  of  RDG  Model  Constants. 
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Figure  3.7.  Comparison  between  Model  Prediction  and  the  Measured  Stress  History 
(OFHC  Copper). 


conventional  plate  impact  experiment  in  that  two  flyer  plates  are  used,  separated  by  a  small  gap, 
as  shown  schematically  in  Figure  3.8.  The  first  flyer  has  a  lower  shock  impedance  than  the 
second  flyer.  When  the  first  flyer  impacts  the  target  plate,  microvoid  nucleation  and  growth 
occur  in  a  localized  area  of  maximum  tensile  stress  in  the  target.  The  reshocking  by  the  second 
flyer  reverses  the  damage  by  closing  the  microvoids. 

Grove  et  al.  [5]  simulated  an  experiment  (829)  in  which  the  first  flyer  was  a  2  mm  thick 
aluminum  plate  and  the  second  flyer  was  a  3  mm  thick  copper  plate.  The  spacing  between  the 
flyers  was  0.25  mm,  and  the  impact  velocity  was  approximately  336  m/s.  The  target  was  a  4  mm 
thick  OFHC  copper  plate  with  a  PMMA  backing.  The  stress  history  at  the  interface  of  the  target 
and  PMMA  was  measured  with  a  manganin  stress  gauge. 

The  double  flyer  plate  impact  experiment  was  simulated  using  the  EPIC-2  finite  element 
code.  EPIC-2  material  models  were  used  to  describe  both  flyers  and  the  PMMA.  The  PMMA 
was  modeled  as  a  rate  independent  elastic-perfectly  plastic  material  with  a  yield  strength  of  0.18 
GPa.  The  RDG  model  was  used  to  describe  the  plastic  flow  and  dilatation  in  the  OFHC  copper 
target  plate.  The  previously  determined  RDG  model  parameters  for  OFHC  copper  were  used  in 
this  simulation. 

In  the  first  simulation  of  the  double  flyer  plate  experiment,  the  RDG  model  predicted  the 
initial  spall  signal  reasonably  well  but  was  unable  to  accurately  reproduce  the  experimental 
recompaction  signal,  as  indicated  by  Figure  3.9.  The  plastic  wave  slope  of  the  simulated 
reshocking  exhibited  ramping  with  a  long  rise  time  indicating  that  the  predicted  pore  collapse, 
described  by  the  same  equations  as  void  growth,  was  unrealistically  slow  and  required 
modification.  Grove  et  al.  [5j  introduced  a  complete  collapse  of  the  microvoids  in  a  particular 
element  if  that  element’s  dilatation  rate  was  negative  and  its  void  volume  fraction  was  less  than 
or  equal  to  some  critical  value.  Figure  3.10  shows  the  excellent  fit  to  the  experiment  obtained 
when  the  critical  void  volume  fraction  for  complete  collapse  was  five  percent. 

The  double  flyer  plate  experiment  without  void  growth  was  also  simulated  so  that  the 
microvoids  effect  on  the  reshock  signal  could  be  investigated.  Figure  3. 1 1  compares  the  two 
simulations,  with  and  without  void  growth.  It  is  interesting  to  note  that  the  dip  at  point  F  occurs 
at  the  same  time  in  both  curves,  indicating  that  this  feature  is  directly  related  to  the  experimental 
configuration.  Figure  3.12  shows  an  x-t  (distance-time)  diagram  for  the  double  flyer  plate  impact 
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Figure  3.8.  Experiment  Configuration  for  Double  Flyer  Impact  Experiment. 
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Figure  3.9.  Stress  History  Comparison  for  Original  RDG  Model. 
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Figure  3.10.  Stress  History  Comparison  for  Modified  RDG  Model,  with  a  Critical  Void 
Volume  Fraction  of  5%  for  Complete  Collapse. 
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Comparison  of  Simulated  Stress  Histories  from  Double  Flyer  Impact 
Configuration,  with  and  without  Void  Growth. 
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Figure  3.12.  X-T  Diagram  for  the  Double  Flyer  Impact  Experiment,  Assuming  No  Void 
Growth. 
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experiment,  assuming  no  void  growth.  In  this  figure,  the  solid  and  dashed  lines  represent  shock 
and  release  waves,  respectively.  The  release  and  shock  waves  that  arrive  at  points  E  and  F  are 
solely  due  to  the  impedance  mismatch  between  the  first  flyer  and  the  target.  Because  the  first 
flyer  has  a  lower  impedance  than  the  target,  it  separates  from  the  target  at  about  0.75 
microseconds.  As  Figure  3.12  indicates,  the  initial  shock  wave  from  the  second  flyer  reflects  off 
the  free  surface  of  the  first  flyer  as  a  release  wave  because  of  the  gap  formed  between  the  first 
flyer  and  the  target.  This  release  wave  can  be  traced  to  point  E.  When  the  first  flyer  impacts 
the  target  again,  another  shock  wave  is  produced,  and  this  shock  wave  can  be  traced  to  point  F 
in  the  x-t  diagram. 

In  the  absence  of  microvoid  nucleation  and  growth,  the  shock  wave  from  the  second  flyer 
would  arrive  at  point  C.  When  void  growth  occurs,  however,  this  signal  is  delayed  because  of 
the  reduced  wave  speed  in  the  porous  region  of  the  target.  Arrival  of  this  signal  at  point  C /  (in 
Figure  3.11)  signifies  that  the  microvoids  have  completely  collapsed.  Since  the  release  and  shock 
waves  at  points  E  and  F  arrive  after  point  C 1 ,  it  appears  that  they  are  unaffected  by  the  void 
growth  phenomenon. 

In  summary,  a  double  flyer  plate  impact  experiment  was  successfully  simulated  using  the 
RDG  model.  Using  previously  calibrated  constants,  the  RDG  model  was  able  to  reproduce  the 
initial  spall  signal,  but  the  reshock  signal  did  not  compare  as  well  with  the  experimental  stress 
history.  The  RDG  model  was  then  modified  to  include  a  simple  void  collapse  criterion  in  which 
the  microvoids  were  forced  to  completely  collapse  below  a  critical  void  volume  fraction  when 
the  dilatation  rate  was  negative.  Microvoid  collapse  in  OFHC  copper  was  successfully  modeled 
using  precalibrated  RDG  model  constants  and  a  critical  void  volume  fraction  for  complete 
collapse  of  5  percent.  Thus,  the  RDG  model  has  been  extended  to  describe  microvoid  collapse 
as  well  as  a  history  dependent  failure  process. 

3.6  RDG  Model  Applications 

Dynamic  ductile  failure  also  occurs  in  other  types  of  configurations.  Tensile  necking  in 
a  split  Hopkinson  bar  test  specimen,  spallation  in  a  target  plate  impacted  by  a  rod  penetrator,  and 
spallation  in  a  solid  cone  impacted  by  a  thin  plate  are  a  few  examples.  These  configurations  are 
investigated  in  the  following  subsections. 
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3.6.1  Tensile  Necking 

To  further  demonstrate  the  capability  of  the  RDG  model,  several  two-dimensional  axi- 
symmetric  configurations  were  considered.  The  first  case  considered  the  problem  of  dynamic 
tensile  necking.  The  main  objective  of  this  two-dimensional  analysis  was  to  examine  the  stress 
states,  evolution  of  necking,  void  growth  as  influenced  by  the  plastic  deformation,  and  ductility. 
To  initiate  necking,  a  shallow-notched  copper  specimen  of  the  geometry  described  in  Figure  3.13 
was  considered.  The  length  and  uniform  diameter  of  the  specimen  are  8.9  mm  and  3.2  mm, 
respectively.  The  minimum  diameter  at  the  notch  is  2.7  mm.  The  length  represents  the  actual 
gauge  length  of  a  standard  SHB  tensile  specimen. 

The  shallow  notched  specimen  was  modeled  using  the  quadrilateral  element  option  in  the 
EPIC-2  code.  Due  to  specimen  symmetry,  it  is  sufficient  to  model  only  one  quarter  of  the 
specimen  as  shown  in  Figure  3.14.  The  Bodner-Partom  (BP)  model  was  used  to  describe  the 
high  strain  rate  behavior  of  the  copper  specimen.  The  BP  model  constants  for  copper  are  given 
in  Appendix  B.  A  typical  experimentally  measured  SHB  velocity  history  was  applied  to  the 
specimen  in  the  simulation.  The  velocity  was  linearly  increased  from  zero  to  50.8  m/sec  (2000 
inches/sec)  for  40  microseconds,  and  afterwards,  it  was  kept  constant. 

The  RDG  model  was  used  to  describe  the  necking  evolution.  A  mean  effective  plastic 
strain  based  Gaussian  distribution,  as  described  by  Equation  (8),  was  used  for  void  nucleation. 
Since  the  stress  level  under  a  uniaxial  stress  state  is  considerably  lower  than  the  level  under 
uniaxial  strain  conditions  (plate  impact  test),  the  void  nucleation  was  assumed  to  occur  entirely 
due  to  large  plastic  flow  around  the  inclusions  and  oxide  particles.  In  a  uniaxial  tensile  test,  the 
necking  often  progresses  under  significant  plastic  flow  and  a  moderate  triaxial  stress  state. 

In  the  analysis,  the  time  histories  of  effective  stress,  effective  plastic  strain,  void  volume 
fraction,  and  triaxiality  of  the  stress  state  (ratio  of  pressure  to  effective  stress)  were  evaluated  for 
elements  82,  182,  282,  and  382  which  are  shown  in  Figure  3.14.  Variation  of  pressure,  stress 
components,  strength,  and  other  variables  along  a  radius  for  different  time  intervals  have  also 
been  considered  in  the  analysis. 
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Figure  3.13.  Shallow-Notched  Tensile  Specimen  Geometry. 
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Figure  3.14.  Finite  Element  Mesh  for  Shallow-Notched  SHB  Specimen.  The  Velocity  Loading 
History  is  Shown  in  the  Inset. 
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The  stress  triaxiality  parameter  (P/Y)  contours  are  plotted  in  Figure  3.15.  While  this 
parameter  is  increasing  at  the  specimen  center,  the  stress  state  is  uniaxial  (P/Y  is  around  1/3)  in 
the  uniform  region.  The  triaxiality  ratio  becomes  negative  (compressive)  within  region  A  as  can 
be  seen  from  Figure  3.15b.  The  unloading  of  the  uniform  section  (region  B)  and  the  compressive 
loading  at  A  make  the  interpretation  of  the  triaxiality  contours  complicated.  In  the  presence  of 
increased  triaxiality,  the  void  growth  also  increases  as  shown  in  Figure  3.16.  The  void  volume 
fraction  contours  at  t=0,  20,  50,  and  60  microseconds  are  shown.  Initially,  a  small  amount  of 
voids  nucleate  near  the  stress-free  notch  surface.  The  void  distribution  spreads  toward  the 
specimen  center  with  increasing  necking.  By  about  50  microseconds,  the  void  content  at  the 
center  region  increases  rapidly  to  3  to  4  percent.  Due  to  increasing  triaxial  tensile  state,  rapid 
void  growth  occurs  during  the  next  10  microseconds;  by  t=60,  the  necked  region  is  filled  with 
voids,  over  40  percent  located  at  the  center  (see  Figure  3.16d).  Although  the  voids  nucleate  near 
the  surface,  the  failure  initiation  eventually  occurs  at  the  center  of  the  specimen.  Hancock  and 
Brown  [36]  provided  the  physical  evidence  to  support  failure  initiation  at  the  center  in  notched 
tensile  specimens  for  various  notch  geometries  under  quasi-static  loading  conditions. 

The  time  history  plots  of  effective  plastic  strain  and  void  volume  fraction  for  regions  near 
the  neck  (local,  element  382)  and  away  from  the  neck  (uniform,  element  82)  are  given  in  Figure 
3.17.  Initially,  the  plastic  strains  are  identical  at  both  regions.  However,  the  local  strain  at  the 
center  of  the  specimen  starts  deviating  from  the  uniform  strain  slowly  due  to  enhanced  void 
growth.  The  void  nucleation  process  occurs  during  the  first  40  microseconds,  and  rapid  growth 
initiates  around  t=50.  Beyond  this  time,  while  the  local  strain  continues  to  increase  due  to  void 
growth,  the  uniform  strain  at  element  82  (see  Figure  3.17)  reaches  a  maximum  value  of  20 
percent  and  stops  increasing.  This  is  due  to  unloading  of  the  sections  away  from  the  local 
necking  region. 

The  stress-time  histories  for  the  local  and  uniform  regions  are  shown  in  Figure  3.18. 
Also  shown  in  this  figure  is  the  void  volume  fraction  history  for  elements  82  and  382.  Since 
these  two  elements  are  very  close  to  the  axis  of  symmetry,  the  hoop  and  radial  stresses  are 
similar  (Figure  3.18).  For  the  unifo.M  element  (82),  the  hoop  and  radial  stresses  are  zero  as  one 
would  expect.  However,  the  stress  state  in  the  local  element  382  begins  almost  uniaxial  and 
becomes  triaxial  due  to  necking.  The  hoop  and  radial  stresses  increase  as  the  void  growth  occurs 
in  the  local  region  generating  a  high  local  triaxial  stress  state.  Between  40  and  50  microseconds. 
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Figure  3.15.  Stress  Triaxiality  (PA')  at  (a)  t=40  and  (b)  t=60  microseconds. 
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amage  (Void  Volume  Fraction)  Contours  at  (a)  t=0,  (b)  t=20,  (c)  t-50, 
60  microseconds. 
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Figure  3.17.  Time  Histories  of  Effective  Plastic  Strain  and  Void  Volume  Fraction  for  Elements 
82  (Uniform  Section)  and  382  (Local  Section). 


Void  Volume  Fraction,  f 
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Figure  3.18.  Time  History  Plots  for  the  Stress  Components  and  Void  Volume  Fraction. 


this  stress  state  accelerates  the  necking  process.  As  rapid  void  growth  occurs  at  50  microseconds, 
the  axial  as  well  as  hoop  and  radial  stresses  relax  (at  point  A 1 )  in  local  element  382.  The  word 
’relax’  indicates  that  the  stresses  decrease  with  increasing  plastic  strains.  However,  the  axial 
stress  actually  unloads  in  uniform  element  82  (at  point  A  in  Figure  3. 18)  due  to  rapid  localization 
of  the  deformation.  Since  no  void  growth  occurs  in  this  uniform  region,  the  decreasing  stresses 
are  only  due  to  unloading.  In  regions  away  from  the  neck,  the  material  is  intact  without  any 
voids,  and,  therefore,  the  possibility  of  any  strength  degradation  leading  to  stress  relaxation  does 
not  exist.  To  further  demonstrate  these  aspects,  the  effective  stress  versus  effective  plastic  strain 
behavior  for  the  two  locations  is  shown  in  Figure  3.19.  The  void  volume  fraction  history  has 
also  been  included.  The  stress  (strength  degradation)  relaxation  beyond  point  A  can  be  seen  in 
this  figure.  The  plastic  strain  continually  increases  because  of  the  void  growth  induced  plastic 
flow.  As  the  void  volume  fraction  increases,  the  strength  continues  to  decrease. 

The  mean  stress  (pressure)  histories  at  various  locations  are  plotted  in  Figure  3.20.  As 
mentioned  earlier,  82  represents  the  uniform  region  and  382  represents  the  local  region. 
Elements  282  and  182  represent  the  intermediate  regions  (see  Figure  3.14).  The  stress  triaxiality 
(defined  as  the  ratio  of  mean  stress  and  effective  stress)  due  to  the  initial  notch  leads  to  a  higher 
mean  stress  level  in  the  local  region  before  the  occurrence  of  necking.  When  necking  initiates 
around  40  microseconds,  the  regions  (82  and  182)  away  from  the  neck  start  unloading.  While 
the  pressure  decreases  because  of  unloading  in  these  regions,  the  mean  stress  also  relaxes  due 
to  void  growth  in  the  local  regions. 

The  stress  triaxiality  is  plotted  with  respect  to  time  for  elements  82  and  382  in  Figure 

3.21.  The  uniform  region  experiences  behavior  close  to  a  uniaxial  stress  state.  The  ratio  of  mean 
stress  to  effective  stress  was  almost  1/3,  as  one  would  expect.  The  slight  deviation  from  this 
value  before  20  microseconds  is  numerical,  and  later  this  ratio  remains  close  to  1/3.  For  the 
element  in  the  local  region,  the  stress  triaxiality  is  initially  greater  than  1/3  due  to  the  notch. 
However,  as  necking  occurs  between  40  and  50  microseconds,  the  local  triaxiality  rapidly 
increases  as  shown  in  Figure  3.21.  Beyond  point  B ,  rapid  void  growth  leads  to  strength 
degradation  and  failure  initiation.  To  further  understand  the  loading  conditions,  the  loading  paths 
are  shown  at  different  locations  with  respect  to  the  pressure  dependent  yield  surfaces  in  Figure 

3.22.  The  loading  path  for  each  element  (82,  182,  282,  and  382)  is  described  by  plotting  the 
ratio  of  effective  stresses  of  the  aggregate  and  the  fully  dense  matrix  (creff/Ym)  versus  the  ratio 
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Response  for  Elements  82  and  382,  along  with  the  Variation  of  f 
to  Effective  Plastic  Strain.  Note  that  f=0  for  Element  82. 
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it  of  Mean  Stress  for  Elements  82,  182,  282  and  382. 
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Measure  of  Triaxiality  with  Respect  to  Time  for  the  Uniform  (82)  and  Local  (382) 
Sections.  The  Dashed  Line  Represents  the  Uniaxial  Stress  State. 


003 


aths  with  Respect  to  Pressure  Dependent  Yield  Surfaces  in  Elements  82, 
and  382. 


of  aggregate  mean  stress  and  matrix  effective  stress  (3P/Ym).  Since  the  material  is  initially  void- 
free,  the  plastic  yielding  is  pressure  independent  and  the  loading  path  starts  at  aef|/Ym  =  1,  as 
shown  in  Figure  3.22.  The  dotted  lines  represent  the  elliptical  yield  surfaces  for  various  levels 
of  porosity.  In  element  82,  a  limited  amount  of  voids  nucleate  initially  and  grow  very  little.  The 
corresponding  loading  path  is  shown  between  points/l  and  B.  The  loading  remains  on  the  yield 
surface  for  void  volume  fractions  less  than  0.3  percent.  In  a  similar  manner,  the  loading  path 
for  element  182,  which  is  away  from  the  notch  zone,  is  confined  to  a  yield  surface  corresponding 
to  f  =  0.005.  Because  of  enhanced  void  growth  in  the  local  regions,  a  significant  amount  of  void 
growth  occurred,  especially  in  element  382  as  between  points  C  and  D  in  Figure  3.22.  As  the 
loading  path  proceeds  along  the  shrinking  yield  surfaces,  the  material  strength  degrades  to  less 
than  80  to  90  percent  of  the  original  intact  material,  indicating  failure  of  the  material. 

To  complete  the  discussion  of  the  RDG  model  capabilities  as  applied  to  the  dynamic 
necking  problem,  additional  analyses  in  terms  of  snapshots  and  additional  contour  plots  are 
presented.  The  snapshots  are  plots  that  describe  the  distribution  of  a  variable  with  respect  to 
position  (for  example,  radial  distance  from  the  axis-of-symmetry).  In  Figure  3.23,  void  volume 
fraction  snapshots  at  three  different  z  positions  (see  Figure  3.14)  are  plotted.  The  z  position 
describes  the  axial  distance  (along  the  specimen  length)  from  the  minimum  cross  section  (c/s) 
of  the  shallow  notched  tensile  specimen.  Therefore,  z=0  represents  the  minimum  cross  section. 
Since  the  triaxiality  is  largest  at  the  center  of  the  specimen  (r=0  and  z=0),  the  void  volume 
fraction  is  maximum  at  the  center  and  decreases  toward  the  free  surface.  The  present  analysis 
is  for  a  dynamically  stretched  tensile  specimen  with  a  loading  duration  on  the  order  of 
microseconds.  Interestingly,  these  numerical  results  are  very  similar  to  those  of  a  statically 
deformed  specimen.  These  analyses  compared  qualitatively  well  with  the  results  reported  by 
Norris  et  al.  [37].  They  analyzed  quasi-static  necking  deformation  using  a  dynamic  finite 
difference  computer  program.  However,  in  the  present  work,  the  problem  is  considered  to  be 
inherently  dynamic  and  has  been  treated  as  a  shock  wave  propagation  problem.  Due  to  several 
reverberations  of  the  shock  waves,  the  deformation  is  homogeneous,  and  the  problem  is 
dominated  by  inertial  loading  rather  than  the  shock  wave  effects. 

Snapshots  of  normalized  effective  stress  in  the  voided  aggregate  are  given  in  Figure  3.24. 
At  t=20  microseconds,  the  effective  stresses  of  the  matrix  and  the  aggregate  are  almost  the  same 
uue  to  insignificant  void  growth;  therefore,  the  stress  ratio  (ocfl/Ym)  is  one.  When  the  void 
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lial  Distribution  of  Void  Volume  Fraction  at  Three  Different  Positions  along 
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Figure  3.24.  Normalized  Effective  Stress  (Gc[[/Ym)  Distribution  along  the  Radius  for  Different 
z  Positions  and  Times. 


growth  occurs,  the  aggregate  strength  degrades  and  the  stress  ratio  becomes  less  than  one.  By 
40  microseconds,  the  c/s  at  z=0.84  mm  and  the  minimum  c/s  have  developed  significant  voids. 
The  z=0  c/s  degrades  faster  than  the  section  at  0.84  mm.  Since  the  void  volume  fraction  is 
initially  evenly  distributed  across  the  c/s  (not  shown  in  the  figure),  the  reduct>cns  in  strength  are 
almost  uniform.  Between  40  and  60  microseconds,  non-uniform  void  growth  in  the  intermediate 
section  leads  to  rapid  strength  degradation  near  the  axis  (r=0),  as  shown  between  A  and  B  in 
Figure  3.24.  However,  because  of  a  lack  of  stress  triaxiality,  the  void  growth  near  the  surface 
(for  the  z=0.84  c/s)  is  not  significant;  correspondingly,  the  strength  near  the  surface  does  not 
significantly  degrade.  By  60  microseconds,  due  to  enhanced  void  growth  and  void  coalescence, 
the  aggregate  strength  at  the  minimum  c/s  (z=0)  has  decreased  by  more  than  90  percent.  This 
indicates  failure  initiation  leading  to  total  separation  of  the  material. 

For  comparison,  the  necking  process  was  also  simulated  without  voids;  this  simulation 
will  be  referred  to  as  case  2  below.  As  before,  the  Bodner-Partom  model  was  used  to  describe 
the  material  behavior,  and  the  solution  was  conducted  for  60  microseconds.  Since  the  boundary 
conditions  for  the  two  simulations  were  the  same,  the  top  section  was  moved  (pulled)  identically. 
Therefore,  at  any  given  instant,  the  specimen  length  was  the  same  for  both  cases.  Figure  3.25 
compares  the  following  strains  from  the  two  simulations:  (1)  ratio  of  current  length  and  original 
length,  (2)  uniform  effective  plastic  strain  in  element  82  for  the  void-free  material  (case  2),  and 
(3)  uniform  effective  plastic  strain  in  element  82  for  the  void-containing  material  (case  1).  As 
expected,  the  average  strain  based  on  the  overall  specimen  length  was  larger  than  the  uniform 
strains  in  element  82  for  the  two  cases.  While  in  case  2  the  uniform  section  continues  to  deform 
at  60  microseconds,  deformation  no  longer  occurs  in  case  1  due  to  the  void  growth  controlled 
necking  process.  For  this  reason,  the  final  strains  in  the  uniform  region  were  significantly 
reduced  in  the  presence  of  void  growth. 

Plots  of  effective  stress  versus  effective  plastic  strain  in  the  local  element  (382)  are 
shown  in  Figure  3.26  for  both  cases  (with  and  without  voids).  In  the  void-free  case,  since 
material  degradation  does  not  occur,  the  effusive  stress  (strength)  follows  the  true  stress-strain 
curves  generated  by  die  BP  model  for  OFHC  copper.  In  the  void-growth  influenced  necking 
case,  the  decreasing  portion  of  the  stress-strain  curve  (between  points  A  and  B)  is  due  to  strength 
degradation.  Even  though  the  loading  (increasing  plastic  strain)  increases  beyond  point  A,  the 
stress  continually  drops,  indicating  void-softening  of  the  material. 


3-41 


NlVtfiS  DliSVld  3AIXD3dda  'NIYdXS  DAV 


3-42 


Figure  3.25.  Comparison  of  Strains  in  Element  82  for  Tensile  Necking  with  and  without  Voids. 

The  Solid  Line  Without  Symbols  is  an  Average  Measure  of  Strain  Based  on  the 
Pull  Velocity. 
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The  effective  plastic  strain  variations  along  a  radius  for  the  necking  process  with  and 
without  voids  are  shown  in  Figure  3.27.  The  variation  along  the  minimum  radius  (local  region) 
indicates  that  the  strain  reached  a  maximum  near  the  surface  and  not  at  the  axis  for  both  cases. 
These  results  differ  from  the  tensile  simulation  results  reported  by  Norris  et  al.  [37]  in  which  the 
plastic  strain  reached  a  maximum  at  the  axis.  They  considered  a  quasi-static  necking  problem 
in  which  the  inertial  effects  were  neglected.  Since  the  present  necking  analysis  considers  the 
inertial  effects,  the  difference  in  the  effective  plastic  strain  distribution  is  attributed  to  the  inertial 
loading.  Also,  the  stress-strain  response  was  modeled  as  an  isothermal  process.  In  the 
simulations,  the  effects  due  to  the  adiabatic  process  are  not  considered. 

The  local  plastic  strain  levels  in  the  void-growth  case  are  larger  than  those  in  the  void- 
free  case.  The  obvious  reason  is  that  void  growth  and  pressure-dependent  yielding  enhance  the 
plastic  flow.  While  loading  continues  in  the  local  region,  the  uniform  section  unloads  because 
of  the  enhanced  necking  process  due  to  void  growth.  In  the  void-free  case,  since  the  necking  is 
not  very  severe,  the  loading  continues  in  the  uniform  section,  as  shown  earlier  in  Figure  3.26. 
Because  of  this,  the  effective  plastic  strain  also  continues  to  increase  in  the  uniform  section. 

Snapshots  of  stress  triaxiality  for  the  two  cases  are  plotted  in  Figure  3.28.  In  the  void 
growth  case,  the  triaxiality  was  maximum  at  the  local  (minimum  radius)  section,  especially  near 
the  axis  (r=0)  as  can  be  seen  from  the  figure.  Other  investigators  [37-39]  have  observed  similar 
behavior  under  quasi-static  loading  conditions,  with  and  without  void  growth.  The  triaxiality 
levels  are  relatively  low  near  the  surface,  as  expected.  In  the  void-free  case,  the  stress  triaxiality 
is  equal  to  about  0.3  in  the  uniform  section.  However,  in  the  uniform  section  of  the  void-growth 
case,  the  triaxiality  is  actually  lower  than  0.3.  The  reason  for  this  can  be  explained  by  Figure 
3.29.  This  figure  shows  a  mean  stress  contour  plot  for  the  void-growth  case.  The  material  under 
compression  is  shown  by  the  region  which  is  not  shaded.  The  contour  was  plotted  for  the  final 
solution  time  (60  microseconds).  Failure  initiation  has  occurred  due  to  a  large  void  content  near 
the  central  portion  of  the  tensile  specimen.  This  final  phase  in  the  failure  process  results  in 
unloading  of  the  uniform  section  and  further  loading  under  compression.  In  the  void-free  case 
(not  shown  in  the  figure),  the  uniform  section  unloads  completely  but  does  not  go  into 
compression.  If  calculations  were  performed  beyond  60  microseconds,  compressive  loading 
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Figure  3.27.  Radial  Distributions  of  e[?ff  for  Cases  with  and  without  Voids  at  t=60  ps,  Plotted 
for  Uniform  (z=1.61)  and  Local  (z=0.0)  Sections. 
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Figure  3.29.  Pressure  Contour  Plot  for  Void-Growth  Influenced  Necking  Process 
The  Shaded  Region  is  under  Tensile  Loading,  and  the  Rest  of 
Experiences  Compressive  Loading. 


might  also  occur  in  the  void-free  case.  However,  Norris  et  al.  [37],  in  their  numerical  analysis 
of  quasi-static  tensile  necking,  observed  compressive  zones  away  from  the  local  region. 

The  deformed  configurations  of  the  shallow  notched  tensile  specimen  for  the  cases  with 
and  without  void  growth  are  shown  in  Figure  3.30.  Only  one  half  of  the  specimen  is  shown  for 
each  case  for  comparison.  Also  included  is  the  initial  geometry,  shown  by  the  dotted  line.  In 
the  void-free  case,  the  neck  profile  retained  the  original  notch  shape,  while  in  the  void-growth 
case,  the  profile  resulted  in  a  sharp  localized  notch  at  the  z=0  c/s. 

In  the  present  application,  the  material  was  modeled  as  strain  rate  dependent  and  highly 
strain  hardening.  Additional  exercises  with  different  material  behaviors  such  as  perfectly  plastic, 
strain-rate  independent,  and  thermal  softening,  with  and  without  voids,  will  shed  additional  light 
on  the  evolution  of  necking. 

3.6.2  Rod  Penetration 

The  main  objective  of  this  simulation  was  to  show  the  spall  evolution  in  a  target  using 
the  RDG  ductile  failure  model.  As  shown  above,  the  model  successfully  describes  the  spallation 
phenomenon  in  the  plate  impact  configuration  and  the  behavior  of  void-growth  enhanced  tensile 
necking  problems. 

An  EPIC-2  code  simulation  of  a  2  inch  long,  1  inch  diameter  solid  copper  rod  impacting 
a  6  inch  diameter,  1  inch  thick  HY100  steel  plate  target  at  an  impact  velocity  of  8333  ft/sec  was 
considered.  The  corresponding  finite  element  mesh  is  shown  in  Figure  3.31.  Material  behavior 
was  modeled  using  the  Johnson-Cook  (JC)  constitutive  model  for  the  copper  rod  and  the  RDG 
model  (with  the  Bodner-Partom  (BP)  constitutive  model)  for  the  HY100  steel  plate.  The 
constants  for  the  BP,  JC,  and  RDG  models  are  given  in  Appendix  B  and  in  Table  3.1. 

The  numerical  solution  was  obtained  for  20  microseconds.  For  analysis,  an  element 
(736)  near  the  axis  of  symmetry  of  the  target  and  at  a  distance  of  about  0.8  inches  from  the 
impact  point  was  examined.  Time  histories  of  pressure,  stress  components,  effective  stress,  and 
damage  (void  volume  fraction)  in  this  element  were  stored  for  post-processing.  Snapshot  plots 
were  also  used  to  display  the  distributions  of  these  variables  along  the  target  radius  at  different 
depths  for  various  times.  Finally,  damage  contours  were  plotted  at  different  times  to  examine 
the  evolution  of  a  spall  zone  inside  the  target. 
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Figure  3.31.  Finite  Element  Mesli  for  the  Rod  Penetration  Problem. 


In  Figure  3.32,  the  normalized  pressure  (mean  stress/aggregate  strength)  is  plotted  with 
respect  to  time  for  element  736  in  the  target.  The  evolution  of  damage  in  terms  of  void  volume 
fraction  is  also  shown.  The  normalized  pressure  is  a  measure  of  the  triaxiality  of  the  stress  state. 
A  negative  normalized  pressure  indicates  that  the  material  is  under  compression,  and  void  growth 
cannot  occur.  To  generate  voids  in  an  intact  ductile  metal,  a  high  triaxial  stress  state  under 
tension  is  required.  The  derivations  of  Rice  and  Tracey  [39]  and  McClintock  [40],  and  the  quasi¬ 
static  experiments  of  Hancock  and  MacKenzie  [41]  support  the  concept  that  a  high  triaxial  tensile 
state  enhances  void  nucleation  and  growth.  The  RDG  model,  which  is  based  on  this  concept, 
was  able  to  model  the  void  growth  accordingly  as  can  be  seen  from  Figure  3.32.  In  a  tensile  test, 
before  the  onset  of  necking,  the  triaxiality  ratio  is  1/3.  However,  in  a  penetration  experiment 
configuration,  the  triaxiality  can  become  large  enough  to  generate  spallation  in  the  target.  In  this 
particular  simulation,  the  tensile  triaxiality  became  extremely  high  (>10  in  Figure  3.32). 

In  element  736,  the  release  wave  from  the  edge  (point  A  in  Figure  3.31)  arrives  at  about 
4  microseconds  and  releases  the  compressive  pressure.  Note  that  the  distance  of  this  element 
from  the  edge  A  in  Figure  3.31  is  around  0.9  inches,  and  the  wave  speed  in  HY100  steel  is  about 
0.23  inches/microsecond.  Then,  at  about  5  microseconds,  the  release  wave  from  the  stress  free 
back  surface  of  the  target  reaches  element  736.  (Note  that  the  initial  shock  wave  will  reflect  as 
a  tensile  wave.)  This  tensile  wave  puts  the  element  under  a  high  triaxial  stress  state  leading  to 
void  nucleation  and  growth.  The  void  volume  fraction  increases  while  the  tensile  stress  state 
prevails.  As  the  material  degrades  due  to  voids,  the  pressure  as  well  as  the  strength  relaxes,  as 
shown  in  Figure  3.32. 

The  strengths  of  the  matrix  material  (intact  HY100  steel)  and  the  porous  aggregate 
(HY100  steel  with  porosity)  for  element  736  are  shown  in  Figure  3.33  with  respect  to  time.  For 
further  understanding  and  evaluation  of  the  RDG  model  capabilities,  the  evolution  of  void  volume 
fraction  (damage)  has  also  been  included  in  the  plot.  Initially,  the  material  is  intact;  therefore, 
the  aggregate  and  matrix  strengths  are  the  same,  and  the  curves  between  points  A  and  B  are 
identical.  Void  nucleation  occurs  at  around  t  =  4.9  microseconds.  Subsequent  void  growth 
degrades  the  aggregate  strength,  as  shown  by  the  sudden  drop  of  the  effective  stress  (Ya)  to  point 
C,  and  the  aggregate  strength  remains  significantly  lower  than  the  matrix  strength  (as  shown  by 
the  dotted  line). 
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Figure  3.32.  Time  Histories  of  Triaxiality  (Mean  Stress/Aggregate  Effective  Stress)  and 


Figure  3.33.  Effect  of  Void  Volume  Fraction,  f  (Porosity)  on  the  Aggregate  Strength  in 
Element  736.  The  Matrix  Strength  (Ym)  is  Independent  of  f. 


In  Figure  3.34,  the  loading  path  in  element  736  was  plotted  with  respect  to  the  pressure 
dependent  yield  surfaces.  Yield  surfaces  for  various  values  of  the  void  content  are  indicated  by 
the  dashed  lines  in  this  Figure.  These  surfaces  are  generated  from  the  relationship  (Equation  (1)) 
between  the  normalized  pressure  and  normalized  effective  strength.  The  pressure  and  the 
effective  stress  are  both  normalized  with  respect  to  the  matrix  flow  strength.  Initially,  the 
material  is  void  free  and,  therefore,  the  effective  stresses  (strengths)  of  the  aggregate  and  the 
matrix  are  the  same.  For  this  reason,  the  loading  path  starts  at  point  A,  where  the  value  of  the 
normalized  effective  stress  is  one.  When  small  amounts  of  voids  (f<0.005)  nucleate  between 
points  A  and  B,  the  pressure  dependence  of  the  yielding  makes  the  aggregate  flow  strength  a 
function  of  pressure.  The  void-containing  aggregate  strength  is  lower  than  the  matrix  strength, 
so  the  ratio  between  these  two  strengths  becomes  less  than  one.  As  the  loading  progresses,  rapid 
void  growth  occurs  between  points  B  and  D.  When  the  void  content  in  element  736  has 
increased  to  about  0.02  (at  C),  the  aggregate  flow  strength  has  dropped  to  about  60  percent  of 
the  original  intact  material  (matrix)  strength.  Further  void  growth  between  points  C  and  D 
degrades  the  aggregate  strength  to  about  10  percent  of  its  original  value.  Beyond  point  D ,  the 
material  completely  loses  its  load  carrying  capacity  in  tension,  and  the  flow  strength,  as  well  as 
the  mean  stress  (pressure),  relaxes  to  zero. 

The  RDG  model  does  not  include  a  void  coalescence  model;  however,  the  void  growth 
model  has  been  formulated  in  such  a  manner  that  rapid  void  growth  occurs  at  a  certain  void 
volume  fraction  level,  indicating  void  coalescence.  This  can  be  clearly  seen  in  Figure  3.34 
between  points  D  and  F. 

The  snapshots  in  Figure  3.35  illustrate  the  radial  extent  of  void  growth  and  aggregate 
strength  degradation  inside  the  target  at  a  depth  of  about  0.8  inches  (same  as  element  736).  In 
this  figure,  radial  distributions  of  void  volume  fraction  and  aggregate  strength  are  plotted  for  two 
different  times  (5  and  20  microseconds).  A  small  number  of  voids  have  nucleated  and  grown 
at  the  center  of  the  specimen  between  A  and  B  at  t=5  ps.  The  conesponding  aggregate  strength 
relaxation  can  be  seen  between  A  '  and  B  1 .  Beyond  a  radial  distance  of  0.5  inches,  the  material 
is  void-free.  By  20  ps,  the  radial  extent  of  the  void  distribution  has  .  pread  to  about  1.75  inches 
(between  C  and  D).  This  results  in  further  degradation  of  the  aggregate  strength,  as  indicated 
by  the  solid  line  between  C  1  and  D  1 . 
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The  void  volume  fraction  contours  for  different  time  intervals  are  plotted  in  Figure  3.36. 
At  5  microseconds,  void  nucleation  and  growth  has  occurred  at  a  distance  of  0.9  inches  from  the 
top  of  the  target.  Since  the  material  in  this  region  experiences  relatively  large  tensile  stresses, 
spall  initiates  there.  The  void  growth  rate  is  relatively  lower  outside  the  region  enclosed  by  the 
contour  in  Figure  3.36a.  Between  5  and  10  microseconds,  the  void  growth  extends  to  a  larger 
region  as  shown  in  Figure  3.36b.  At  15  microseconds,  the  spall  process  is  complete  and  intense 
void  concentration  has  occurred.  The  inner  contours  in  Figure  3.36c  represent  this  intense  spall 
zone.  The  void  volume  fraction  in  this  zone  was  well  above  0.5.  In  a  penetration  experiment, 
the  material  often  physically  separates  and  forms  stress  free  fracture  surfaces  in  the  target  due 
to  void  coalescence;  in  some  cases,  spall  fragments  eject  out  from  the  back  of  the  target. 

For  completeness,  the  case  of  the  rod  penetration  process  without  any  material 
degradation  was  simulated.  This  calculation  did  not  include  the  effects  of  damage  or  a  spall 
criterion.  The  material  behavior  was  therefore  described  by  the  JC  or  BP  viscoplastic  models 
and  without  the  RDG  model.  In  Figure  3.37,  the  matrix  and  aggregate  strength  time  histories 
from  the  two  analyses,  one  with  material  degradation  (RDG  failure  model)  and  the  other  without 
any  degradation  model,  are  plotted.  Again,  these  plots  are  for  element  736  of  the  finite  element 
mesh  shown  earlier  in  Figure  3.31.  Since  the  aggregate  and  matrix  strengths  are  the  same  in  the 
absence  of  any  porosity,  the  strength  remained  high  (in  the  case  without  voids)  as  represented 
by  the  dashed  line  in  Figure  3.37. 

The  deformed  configurations  of  the  penetrator  and  the  target  at  time  20  microseconds  are 
compared  in  Figure  3.38.  The  target  exhibited  a  relatively  larger  bulge  from  enhanced  plastic 
flow  due  to  void  softening  of  the  material.  Since  the  damaged  target  material  in  front  of  the 
projectile  allows  the  projectile  to  penetrate  more  easily,  the  penetration  process  is  accelerated. 
When  there  is  no  material  degradation  due  to  void  nucleation  and  growth,  the  target  material  is 
relatively  stronger  and,  therefore,  the  penetration  rate  is  slower  compared  to  the  case  in  which 
the  material  degrades.  The  differences  in  the  deformed  target  shapes  and  the  penetration  rates 
demonstrate  the  importance  of  an  accurate  dynamic  failure  model. 
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Figure  3.37.  Comparison  of  Aggregate  Yield  Strengths  (Ya)  with  and  without  Voids  in  Element 
736.  Ym  is  the  Matrix  Strength  of  the  Void  Containing  Material.  Note  that  Yra 
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Figure  3.38.  Deformed  Configurations  of  the  Rod  Penetration  at  20  p  sec,  with  and  without 
Void  Growth. 


3.6.3  Spall  in  a  Solid  Cone  Target 

Spall  under  a  three-dimensional  strain  state  was  modeled  using  the  RDG  ductile  failure 
model  [2,6].  Using  this  model,  we  simulated  a  plate  impact  configuration  in  which  a  flyer  plate 
impacts  the  base  of  a  solid  right  circular  cone  as  schematically  shown  in  Figure  3.39a.  The  main 
objective  of  this  simulation  was  to  describe  spall  under  a  three-dimensional  strain  state  using  the 
RDG  model  constants  determined  from  a  conventional  plate  impact  experiment. 

We  considered  an  experimental  configuration  which  involved  the  impact  of  a  thin  circular 
plate  on  the  base  of  a  right  circular  cone.  The  flyer  plate  was  2  mm  thick  and  38  mm  in 
diameter,  and  the  cone  height  was  25  mm.  Both  the  flyer  and  target  were  1020  steel.  The 
experimentally  observed  spall  patterns  inside  the  cone  are  shown  in  Figure  3.39b.  It  can  be  seen 
that  four  separate  fracture  systems  developed  in  the  target.  A  one-dimensional  spall-type  fracture 
developed  parallel  to  the  cone  surfaces.  This  fracture  was  apparently  the  first  one  to  develop  and 
it  is  first  manifest  at  about  mid-height.  The  second  fracture  system  consisted  of  distributed 
cavities  near  the  axis.  The  third  system  below  mid  height  toward  the  base  is  roughly  a 
cylindrical  failure  region  about  12  mm  in  diameter  symmetric  to  the  cone  axis.  The  fourth 
system  consisted  of  radial  cracks.  This  failure  is  believed  to  have  occurred  last,  as  a  result  of 
surface  motion  induced  by  wave  reflection  from  the  flanks  of  the  cone. 

The  RDG  model  constants  were  calibrated  by  matching  the  computed  stress  history  with 
the  plate  impact  test  data  as  shown  earlier  in  Figure  3.4  for  1020  steel.  The  cone  impact  problem 
was  simulated  using  the  EPIC-2  code.  The  RDG  model  results  showed  the  evolution  of  the 
experimentally  observed  fracture  regions  (patterns)  as  can  be  seen  in  Figure  3.40.  The  spall 
region  parallel  to  the  cone  surface  evolved  clearly  in  the  simulation.  The  model  also  reproduced 
the  second  region  containing  the  distributed  voids.  However,  the  evolution  and  location  of  the 
cylindrical  spall  zone  was  not  clear. 

For  further  understanding  of  the  code  results,  we  plotted  the  effective  plastic  strain 
contours  in  Figure  3.41.  These  contours  indicate  the  evolution  of  a  cylindrical  spall  plane 
parallel  to  the  cone’s  axis.  To  generate  the  (physical)  spall  planes  in  the  code  simulation, 
elements  were  failed  according  to  a  failure  strain  option  in  the  EPIC-2  code.  Using  this  option, 
elements  that  experienced  a  critical  amount  of  effective  plastic  strain  were  removed  from  the 
finite  element  mesh.  Figure  3.42  shows  the  code  generated  spall  plane  parallel  to  the  lateral 
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surface,  which  compares  well  with  the  experiment  (see  Figure  3.39).  Even  though  the  cylindrical 
spall  plane  did  not  appear  in  Figure  3.42,  the  effective  plasdc  strain  contours  (Figure  3.41) 
confirmed  the  RDG  model’s  ability  to  predict  the  development  of  a  cylindrical  spall  plane.  Using 
a  smaller  critical  failure  strain  to  produce  the  cylindrical  spall  plane  resulted  in  excessive  failure 
at  the  base  o*  fhe  cone. 

As  a  final  exercise,  this  simulation  was  performed  without  void  nucleation  and  growth. 
Figure  3.43  shows  the  resulting  effective  plastic  strain  contours.  Comparing  Figure  3.43  with 
Figure  3.41,  it  is  obvious  that  the  void  nucleation  and  growth  process  has  a  significant  influence 
on  the  evolution  of  plastic  strain  regions  inside  the  cone.  Therefore,  to  successfully  predict  spall, 
it  is  important  to  use  models  that  incorporate  degradation  of  stiffness  and  strength  due  to  void 
nucleation  and  growth. 
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Section  4 

Impact  Damage  Model  for  Ceramics 


Ceramic  materials  exhibit  inelastic  deformation  when  shocked  above  the  Hugoniot  Elastic 
Limit  (HEL).  Due  to  lack  of  any  strong  physical  evidence,  it  can  only  be  speculated  that  the 
inelastic  deformation  is  due  to  either  plasticity  (dislocations  movement)  or  microcracking  or  both. 
The  microscopic  examination  of  particle  impacted  ceramic  targets  often  reveals  dislocation 
controlled  plastic  flow  at  the  impact  site.  Such  evidence  has  not  been  established  yet  in  flyer 
plate  impacted  targets,  especially  at  impact  velocities  above  the  HEL.  However,  in  a  few  ceramic 
materials,  microcracks  have  been  observed  in  targets  (plates)  impacted  below  the  HEL  level.  In 
1988,  Rajendran  and  Cook  [9]  summarized  a  literature  review  on  the  impact  behavior  of 
ceramics. 

When  developing  the  constitutive  model  to  describe  the  impact  behavior  of  ceramic 
materials,  a  microphysical  based  approach  was  chosen,  instead  of  an  empirical  or 
phenomenological  approach,  because  it  seemed  best  suited  for  the  analysis  of  both  the 
microcracking  and  plastic  processes.  Typically,  a  microphysical  model  formulation  begins  by 
decomposing  the  strain  into  its  elastic  and  inelastic  components.  Physically  motivated  equations 
are  then  used  to  relate  the  inelastic  strains  to  underlying  mechanisms  such  as  microcracking  and 
void  opening/collapse  (plasticity).  Since  ceramics  fail  at  low  stresses  under  tensile  loading, 
inelastic  strains  are  assumed  to  be  strictly  due  to  microcracking  only.  The  inelastic  strains  due 
to  dislocation  motion  can  be  described  by  viscoplastic  equations. 

In  the  ceramic  model  developed  here,  Margolin’s  equations  [12]  for  a  penny  shaped  crack 
are  used  to  calculate  the  inelastic  strains  due  to  crack  opening  and  sliding  both  under  tension  and 
shear  loading.  However,  under  compression,  only  the  displacement  fields  due  to  crack  sliding 
are  considered.  Therefore,  crack  opening  under  compression  is  not  considered  in  the  present 
formulation.  Further,  ceramic  materials  are  assumed  to  contain  a  large  number  of  closed 
microflaws  and  that  the  microcracks  propagate  upon  shock  loading.  The  model  formulation 
under  this  assumption  eliminates  any  nucleation  model  or  threshold  condition  requirements;  this 
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minimizes  the  number  of  model  parameters.  This  section  presents  the  mathematical  formulation 
of  this  microphysical  based  model  for  ceramic  materials. 

4.1  Constitutive  Equations 


In  the  case  of  a  ceramic  containing  microcracks  and  spherical  voids,  the  total  strain  may 
be  decomposed  into  elastic  and  plastic  strains  as, 


e  .  ,  =  ee- +e?. 
ci  j  e  ij 


where  the  elastic  strain  consists  of  the  matrix  elastic  strain,  the  microcrack  opening  strain,  and 
the  elastic  strain  of  the  embedded  spherical  voids, 


e  m  c  v 

p  .  .  s  p  •  ‘  +  p  .  .  +  p  .  . 
*  IJ  fc2J  blJ 


(15) 


The  matrix  elastic  strain  e™j  and  the  microcrack  opening  strain  e  Jj  are  both  proportional  to  the 

applied  stress  field.  By  definition,  the  total  elastic  strains  of  the  spherical  voids  e  ^  can  be 
expressed  as. 


„  V  _  v  .  ev«? 
eij  ~  eij  y  °ij 


(16) 


where  ejj  are  the  elastic  deviatoric  strains  and  e  y  is  the  elastic  volumetric  strain  associated  with 
the  voids.  These  strains  can  be  defined,  using  Mackenzie’s  equations  [30],  as 


oV  —  pe  —  R  Pe 
*•  V  ~  “  V  V 


(17) 


pY.  S  of .  -  R  pe. 
ci  j  cij  ™gcij 


where 


B  -  d-f ) 

Hk  ~  -) - = — 

l+2Kf 
4  G 


(19) 
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and 


Rg=(l~f) 


i  _  (6K  +  1 

9  K  + 


12  G)  f 
+  8  G 


(20) 


K  is  the  bulk  modulus  and  G  is  the  shear  modulus  of  the  microcracked  material  without  spherical 
voids,  and  f  is  the  void  volume  fraction  related  to  the  amount  of  volumetric  plastic  strain. 
Using  Equations  (14)  through  (18),  the  elastic  strains  in  the  cracked  "matrix"  material  may  be 
expressed  as. 


®  v  ^  ®  i;  =  v  “  ^Jr  (  ®  v  ®  v) 


(21) 


eij  +  eij  ~  Rgeij  ~  Rg  (  ei  j  eij )  ' 


(22) 


where  ep  =  e?1  and  efj=  (e?j-ep/ 3)  by  definition. 


A  pressure  dependent  yield  function  to  define  plastic  flow  in  the  ceramic  material  is 
employed.  This  yield  function  has  an  implicit  analytic  form. 


<b(f,  Ym,  ai7)  =0 


(23) 


where  f  is  the  void  volume  fraction,  Ym  is  the  flow  stress  in  the  matrix  (intact)  material,  and 
O;  7-  are  the  total  stresses.  The  total  plastic  strain  rates  are  given  by: 


(1  -f)YmD, 


p  3d> 


tP .  = 


dCTv-7 


where 


f=  1  -  exp  (-®  v) 


(24) 


(25) 


WP  -  Y  DP 


(26) 
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and 


D%  =  g(YmtDP,wP)  . 


(27) 


In  the  above  equations,  Ym,  W%,  and  represent  flow  stress,  plastic  work,  and  equivalent 
plastic  strain,  respectively,  in  the  matrix  (void-free)  material.  The  value  of  Ym  is  obtained  by 
solving  Equation  (23),  using  estimated  values  of  the  void  volume  fraction  and  the  stresses.  The 
plastic  strains  are  obtained  from  the  numerical  integration  of  Equation  (24).  Equation  (25) 

defines  the  relationship  between  the  void  volume  fraction  f  and  the  volumetric  plastic  strain  e 

Equation  (26)  defines  the  plastic  work  rate,  and  Equation  (27)  is  the  flow  rule  for  the  yield 
function;  either  the  Johnson-Cook  or  the  Bodner-Partom  model  can  be  selected.  The  details  of 
the  Equations  (23)  through  (27)  can  be  found  in  Appendix  C. 

To  implement  the  constitutive  model  in  a  hydrodynamic  code,  the  stress  tensor  must  be 
expressed  as  a  direct  function  of  the  total  strain  tensor.  First,  assume  that  the  elastic  strains  in 
the  microcracked  "matrix"  (void-free)  material  are  related  to  the  stress  tensor  as  follows: 


*•  ij  *■  ij 


~  cijkl  akl 


(28) 


where  Cijkl  is  the  effective  compliance  tensor  of  the  microcracked  material.  If  Cijkl  can  be 
analytically  inverted  to  the  stiffness  tensor  Mijkl,  the  resulting  stress  state  is, 

aij  ~  Mijkl  ^ekl+ekl'> 

Assuming  that  the  stiffness  tensor  is  isotropic,  Equation  (29)  can  be  decomposed  by  standard 
procedures  into  pressure  and  deviatoric  stress  expressions: 

P=  -K{z™  +  €y  )  (30) 


and 

Sij  =  2  g" ( efj  +  efj  )  . 

Substituting  Equations  (21)  and  (22)  into  Equations  (30)  and  (31), 

P=  -RkK(ev-e%) 


(31) 


(32) 
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and 


Sjj  -2  Rg  G  (e^j- efj  )  .  (33) 

Recall  that  the  total  plastic  strains  are  available  from  Equation  (24).  The  form  of  Equation  (32) 
suggests  that  the  Hugoniot  equation  of  state  in  the  microcracked  void-free  material  can  be 
approximated  by, 

PH=  (K/ K)  PH  (34) 

where 

^=PiH  +  P242  +P3H3  (35> 


and 


H  =  exp[-(ev-e£)  ]  -1  .  (36) 

Equation  (36)  defines  the  elastic  compressibility  J!  as  a  function  of  the  true  elastic  volumetric 
strain.  The  Mie-Gruneisen  adjustment  for  the  shock  jump  condition  within  the  microcracked 
material  leads  to  the  following  equation  of  state: 

P  =  Rk  [Fh(1-0.5I»  +  rp0(I-I0)]  (37) 

where  T  is  the  Mie-Gruneisen  parameter,  p0  is  the  material’s  initial  density,  IQ  is  the  initial 
value  of  internal  energy,  and  I  is  the  current  internal  energy.  Note  that  Equation  (37) 
incorporates  Mackenzie’s  adjustment  factor  Rk  for  non-zero  values  of  f .  If  T  is  zero  and  the 
elastic  compressibility  p.  is  negligible.  Equation  (37)  reduces  to  Equation  (32).  Similarly,  if  there 
are  no  microcracks  or  voids.  Equation  (37)  reduces  to  the  usual  Mie-Gruneisen  equation  of  state 
for  an  undamaged,  flawless  material.  The  appearance  of  the  internal  energy  term  in  Equation 
(37)  requires  the  numerical  integration  of  the  internal  energy  rate, 

Poi=  (38) 

ro 
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4.2  Damaged  Stiffness  Tensor 

It  is  essential  that  the  damaged  compliance  tensor  be  analytically  invertible  for  practical 
use  of  the  model  in  explicit  hydrodynamic  codes.  This  requirement  of  analytical  formulae  for 
the  damaged  stiffness  tensor  Mijkl  permits  only  isotropic  or  orthotropic  forms  of  Cijkl. 
Because  of  this,  the  present  ceramic  model  formulation  adopts  the  elastic  moduli  of  a  cracked 
body  proposed  by  Margolin  [12]  and  by  Budiansky  and  O’Connell  [42].  Both  of  these 
formulations  lead  to  effective  compliance  tensors  that  are  isotropic. 

For  non-interacting,  penny-shaped  microcracks  of  various  sizes  and  in  random 
orientations,  Margolin  provided  the  following  expression  for  the  isotropic  elastic  moduli: 

Cijkl  ~  cl^ik^jl  +  c2^il^jk  +  c3^ij^kl  (39) 


where 


Ci  =  Bc  + 


1 

TG 


(40) 


C2  =  -Do  + 


1 

TG 


/ 


and 


c2=Aq 


V 

2 (1 +v) G  * 


(41) 


(42) 


In  the  above  equations,  G  and  v  are  the  shear  modulus  and  Poisson’s  ratio,  respectively,  of  the 
undamaged  material,  while  A0,  B0,  and  D0  are  damage  parameters  whose  values  depend  on  the 
stress  state.  To  evaluate  these  parameters,  Margolin  defined  a  microcrack  density  parameter. 


y* -  16 Y 

7  T5B  ' 


(43) 


where 


da  ~  No 


amax  • 


(44) 
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In  Equation  (43),  y  is  the  microcrack  density  and  E  is  the  Young’s  modulus  of  the  undamaged 
material.  Equation  (44)  is  the  expression  for  microcrack  density  that  results  from  assuming  an 
exponential  size  distribution  of  microcracks.  N0  is  the  number  of  microcracks  per  unit  volume,  a" 
is  the  average  microcrack  size,  and  amax  is  the  maximum  microcrack  size  in  the  distribution. 
Margolin  identified  the  following  four  cases  of  stress  state  in  evaluating  the  damage  parameters 
Aol  B0,  and  D0: 


Case  1:  alt  a2,  c3  >  0 


(all  principal  stresses  are  tensile) 

A0  =  ( ( 1 -v2 ) - ( 1+v) ] y* 

(45) 

B0  =  [ ( 1-v2 )  +  4  (1+v)]  Y* 

(46) 

Do  —  AQ 

(47) 

Case  2:  a1 ,  a2,  o3  <  0 

(all  principal  stresses  are  compressive) 

A0  =  - (1+V) y*  (48) 


B0  =  4 (1+v) Y*  (49) 

D0  =  A0  (50) 

Equations  (48)  through  (50)  will  result  in  the  degradation  of  the  shear  modulus,  but  not  of  the 
bulk  modulus,  because  under  compression  only  crack  movement  of  the  closed  microcracks  under 
modes  II  and  III  are  permitted. 
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Case  3:  a1 ,  a2  >  0 ,  a3  <  0 

(two  principal  stresses  are  tensile,  and  one  principal  stress  is  compressive) 


Ap  a  ^(5P3  ~3(35)  (1-y2)  _  (1+V)  y* 


(51) 


3C=  -3(35)  (1  _v2)  +  4(1+V)Jy* 


(52) 


D0  =  [(6p5-5p3)  (1-v2)  -  (1+v)] Y* 


(53) 


where 


B  ss  +®2 

\  c3  +a2  -2cr3 


Equation  (54)  is  an  approximation  to  Margolin’s  numerical  evaluation  of  the  integral  over  solid 
angles  (Equation  6.3  in  Margolin  [12]),  assuming  an  averaged  portion  of  the  solid  angle  is  under 
tensile  stress. 


Case  4:  >  0 ,  a2 ,  a3  <  0 

(one  principal  stress  is  tensile,  and  two  principal  stresses  are  compressive) 
Ap  =  5(1~P3)  ~3(1:Pll(i-v2)  -  (l+v)  jy*  (55) 

R^[5(l-P3)-3(l-P5)  fi_y2i  +  4  (1+V)  ]y*  <56> 


D0  =  {[6(l-P5) -5(1-P3)  ]  (1-v2)  -  (1+v)}  Y*  (57) 
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where 


P 


CT2  +O3 
CT2 


(58) 


As  in  Case  3,  Equation  (58)  is  an  approximation  to  the  numerical  solution  of  the  integral  over 
solid  angles. 

Next,  the  elements  of  the  stiffness  tensor  are  solved  from 


M ijkl  Cijkl  ~ 


(&ik&jl  +  SilSjk) 
3 


(59) 


with  the  result 


Mijkl  ~ml&ikSjl  +m2&il&jk+m3&ijSkl 


(60) 


where 


Ml  =  ™2  = 


1 

2  ( cx  +  c2  ) 


(61) 


and 


77I3  = 


-C3 

( CX  +  c2  )  ( cl  +  c2  +  3  c3  ) 


(62) 


The  degraded  shear  and  bulk  moduli  are  then  given  by 

G  =  mi  (63) 


and 


K-m3 


(64) 
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Budiansky  and  O’Connell  [42]  considered  the  case  where  microcracks  are  randomly 
oriented,  interacting,  self-consistent,  and  subjected  to  tensile  loading.  Relating  their  results  to 
Equation  (60),  the  damaged  stiffness  tensor  solution  is  given  as: 


where 


and 


=  m2  =  G  , 

(65) 

_  2  GV 

(66) 

(1-2V) 

§ 

(67) 

(l-v)  (5-V)  v\ 

(68) 

\  \Tsl 

T^v)  1 J  ' 

In  these  equations,  V  and  G  are  the  Poisson’s  ratio  and  shear  modulus,  respectively,  of  the 
microcrack  damaged  material.  Using  m1  and  m3  from  Equations  (65)  and  (66),  the  degraded 
bulk  modulus  can  be  computed  from  Equation  (64).  It  is  obvious  from  Equations  (65)  through 
(68)  that  a  complete  loss  of  strength  is  predicted  when  the  microcrack  density  y  reaches  9/16. 
For  tensile  loading  conditions,  based  on  comparisons  with  more  detailed  models  (Nemat-Nasser 
and  Obata  [43]),  as  well  as  in  Margolin’s  approach,  there  is  no  bound  on  the  crack  density. 
However,  Budiansky  and  O’Connell’s  solution  limits  the  crack  density  to  9/16.  This  permits  the 
damage  parameter  to  vary  from  zero  (no  damage)  to  one  (fully  damaged).  Therefore,  in  the 
ceramic  model,  Budiansky  and  O’Connell’s  equations  are  used  instead  of  Margolin’s  equations 
for  the  case  when  all  the  principal  stresses  are  positive  (Case  1). 

4.3  Microcrack  Damage  Evolution 

The  microcrack  growth  mechanism  is  analyzed  using  a  dynamic  linear  fracture  mechanics 
theory.  This  mechanism  consists  of:  (1)  nucleating  microcracks  when  the  stress  state  satisfies 
a  generalized  Griffith  criterion,  (2)  propagating  the  cracks  at  a  crack  tip  speed  which  is  a  function 
of  stress  intensity  factor,  and  (3)  modeling  coalescence  of  the  microcracks  at  some  critical  crack 
size.  The  three  stages  of  microcrack  damage  (nucleation,  growth,  and  coalescence)  occur  at 
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different  times  within  a  given  elemental  volume.  Since  current  experimental  diagnostics  are 
incapable  of  resolving  these  stages  of  microcrack  damage,  incorporating  this  process  into  a 
constitutive  theory  becomes  a  challenging  task. 

4.3.1  Microcrack  Nucleation 

In  brittle  materials,  such  as  ceramics,  a  microcrack  starts  to  grow  when  the  stress  field 
at  the  crack  tip  satisfies  Griffith’s  criterion  [44].  If  the  microcracks  in  a  given  material  have  an 
exponential  crack  size  distribution,  then  the  various  microcracks  begin  to  grow  at  different  times 
as  the  stress  levels  increase.  At  low  enough  strain  rates,  the  largest  microcrack  begins  to  grow 
first,  causing  stress  relaxation  in  other  zones,  thus  preventing  other  microcracks  from  growing. 
In  dynamic  plate  impact  situations,  the  strain  rates  are  of  order  103/sec  and  higher.  At  such  high 
strain  rates,  since  the  crack  tip  speed  is  limited  to  the  Rayleigh  velocity,  significant  stress 
relaxation  will  not  occur  until  nearly  all  the  microcracks  have  begun  to  propagate.  In  addition, 
with  the  stress  intensity  factors  becoming  much  higher  than  the  critical  value  in  a  very  short  time, 
multiple  crack  branching  will  occur.  Modeling  the  details  of  these  processes  would  be  an 
extremely  difficult  task. 

Our  simplified  approach  is  the  following.  Microcracks  are  considered  nucleated  (actively 
growing)  when  the  stress  state  is  such  that  Griffith’s  criterion  is  satisfied  for  the  largest 
characteristic  microcrack.  All  smaller  microcracks  are  then  assumed  to  propagate  at  appropriate 
rates  that  maintain  an  exponential  size  distribution  of  microcracks  (in  random  orientation).  These 
assumptions  permit  a  s.  Jghtforward  application  of  linear  fracture  mechanics  concepts  with  some 
caveats.  For  example,  because  it  is  endowed  with  a  statistical  property,  the  characteristic 
microcrack  size  in  the  model  may  not  correspond  to  the  true  microcrack  size  in  the  ceramic. 
Also,  since  the  actual  microcrack  size  distribution  might  not  fit  an  exponential  distribution,  a  bias 
can  occur  in  the  characteristic  crack  size.  If  the  microcracks  are  oriented  in  preferred  directions, 
the  assumption  of  randomly  oriented  cracks  could  introduce  a  bias  to  the  damaged  moduli.  There 
is  the  additional  problem  of  examining  all  crack  orientations  to  determine  the  highest  stress 
intensity  factor  (or  crack  energy  release  rate)  for  nucleating  the  microcracks.  All  these  modeling 
ambiguities  are  tolerable  when  weighed  against  the  requirement  to  develop  analytically  tractable 
solutions  for  microcrack  damage  evolution. 
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To  model  microcrack  nucleation,  a  generalized  Griffith  criterion  developed  by  Margolin 
[12]  and  Dienes  [45]  is  used  for  any  given  crack  orientation.  For  crack  surfaces  perpendicular 
to  tensile  normal  stresses  okk,  the  crack  energy  release  rate  for  mode  I  crack  opening  is 


g+=  4(i,:,y2>  a 

KE 


max 


2  2  (&jl.+a2jk) 

- 


i*j*k 


(69) 


where  amax  is  the  size  of  the  largest  microcrack,  and  v  and  E  are  Poisson’s  ratio  and  Young’s 
modulus,  respectively,  in  the  undamaged  material.  The  repeated  index  ’k’  does  not  mean 
summation.  The  bar  over  the  stress  components  means  stresses  in  the  cracked  material.  When 
the  crack  surfaces  are  perpendicular  to  tensile  principal  stresses  Gj,  then  Equation  (69)  is 
evaluated  with  okk  =  a1  and  Gik  =  Gjk  =  0 .  For  crack  surfaces  perpendicular  to  compressive 
normal  stresses  akk ,  Dienes  derived  the  crack  energy  release  rate  as 


r-_  8 ( 1 -V2 )  _ 

keTT-vT  max 


G7k+Cf2jk  ~xo+(0Okk 


where  x0  is  the  cohesion  stress  and  to  is  the  friction  coefficient.  In  this  case,  only  modes  II  and 
III  are  active  and  the  normal  stresses  serve  only  to  resist  the  shearing  stresses.  Defining  Gmax 

to  be  the  maximum  of  all  values  of  G*  and  G",  microcracks  are  assumed  to  have  nucleated  if 
Gmax  exceeds  the  critical  crack  energy  release  rate  Gc,  where 


Gc  =  2T 


(71) 


and 


T  = 


(72) 


Equation  (72)  defines  the  surface  tension  T  as  a  function  of  the  static  fracture  toughness 
(KIC),  Poisson’s  ratio  (v),  and  Young’s  modulus  (£)  in  the  undamaged  material.  Note  that  all 
three  modes  of  crack  opening  are  considered  in  this  approach. 
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4.3.2  Microcrack  Growth 


From  the  theoiy  of  dynamic  fracture  mechanics  [46],  the  classical  solution  for  the 
dynamic  energy  release  rate  under  mode  I  conditions,  Gj,  is  given  by  the  formula, 


(73) 


where  A  is  the  microcrack  growth  rate.  CR  is  the  Rayleigh  wave  speed,  and  Gc  is  the  critical 
crack  energy  release  rate,  as  defined  by  Equations  (71)  and  (72).  In  generalizing  Equation  (73) 
for  all  three  modes  of  microcrack  growth,  empirical  constants  to  control  the  limiting  crack  tip 
speed  were  introduced.  This  also  allows  a  nonlinear  dependence  on  the  crack  energy  release  rate. 
The  newly  proposed  microcrack  growth  evolution  equation  is  defined  as 
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Qnax  >  Gc 


\  ^max^^c  and  Gmax>Gc 


(74) 


where  G*ax  is  the  maximum  G+  calculated  from  Equation  (69)  and  G^ax  is  the  maximum  G" 
calculated  from  Equation  (70).  Two  distinct  sets  of  microcrack  growth  constants,  one  for  mode 
I  crack  opening  (n*  and  n^)  and  one  for  modes  II  and  III  (n{  and  )  are  considered.  Then* 

and  n{  coefficients  directly  control  the  limiting  crack  growth  rates  for  the  different  modes  of 

crack  opening,  while  the  n^  and  nj  exponents  serve  to  inhibit  crack  growth  during  the  initial 

phase.  Small  positive  values  («1)  of  n^  and  nj  will  effectively  delay  microcrack  growth  until 
the  critical  crack  energy  release  rate  Gc  has  been  significantly  exceeded.  Second  order  effects 
on  the  microcrack  growth,  such  as  toughening  and  crack  interactions,  are  ignored  in  the  model. 
The  microcrack  growth  rate  is  ultimately  a  function  of  the  current  stress  state  and  maximum 
microcrack  size.  Therefore,  the  microcrack  damage  evolution  is  obtained  from  the  numerical 
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integration  of  Equation  (74)  and  the  evaluation  of  the  crack  density  function  from  Equation  (44). 
Since  the  damage  evolution  itself  degrades  the  moduli,  the  stress  levels  will  relax  as  the 
microcracks  grow,  possibly  enough  to  arrest  the  crack  growth.  In  some  cases,  however,  the  stress 
relaxation  will  be  gradual  enough  to  result  in  microcrack  coalescence  before  crack  growth  arrest 
can  occur. 


4.3.3  Microcrack  Coalescence 


Microcrack  coalescence  occurs  when  propagating  microcracks  intersect  each  other  to  the 
point  of  pulverization  and  significant  loss  of  strength.  Based  on  this  definition,  a  reasonable 
pulverization  criterion  can  readily  be  derived  as  follows.  If  N0  is  the  number  of  microcracks  per 
unit  volume,  then  the  average  spacing  between  cracks  is  obtained  by  inverting  Equation  (44),  i.e. 


d-Nj 


(75) 


If  a  is  the  average  crack  radius,  pulverization  can  be  expected  to  occur  when  this  radius  reaches 
half  the  average  spacing,  or 

ap  =  ^  .  (76) 

This  coalescence  condition  defines  a  critical  value  for  the  microcrack  density  function 
as  follows.  Assuming  that  the  constant  ratio  a^ax  is  much  larger  than  unity  and  all  the  flaws 

have  been  transformed  to  penny  shaped  microcracks,  then  a  reevaluation  of  Equation  (44)  results 
in  the  following  approximation  for  y: 

Y  =  K*Lx  “  6N05'3  .  (77) 

Substitution  of  Equations  (75)  and  (76)  into  Equation  (77)  provides  the  following  definition  of 
pulverization  in  terms  of  microcrack  density: 

<78> 


4-14 


This  criterion  is  introduced  into  the  ceramic  model  because  Margolin’s  formulae  for  the  damaged 
moduli  does  not  include  the  effects  of  microcrack  interactions.  The  material  is  considered  to  be 
pulverized  when  Equation  (78)  is  satisfied  (y>  yp)  during  a  stress  state  in  which  at  least  one  of 
the  principal  stresses  is  compressive  (Cases  2,  3,  or  4  from  Section  4.2).  At  this  time,  the 
pulverized  bulk  and  shear  moduli,  and  Gp,  are  defined  to  be  the  current  values  of  K  andG 
computed  from  Equations  (63)  and  (64).  Henceforth,  the  material  has  no  strength  in  tension,  and 
its  compressive  strength  follows  a  Mohr-Coulomb  law,  as  in 


0,  P<  0 

(Xp+PpP,  P>0 


(79) 


where  Y  is  strength,  P  is  pressure,  and  ap  and  (Jp  are  model  constants  for  the  pulverized 
material.  The  pressure  is  computed  simply  from 


0, 


P  = 


-KpC®, 


e*>0 

e  v<  0 


(80) 


where  e  %  is  engineering  elastic  volumetric  strain  and  Kp  is  the  pulverized  bulk  modulus.  With 
this  approach,  each  pulverized  element  in  a  finite  element  mesh  may  have  its  own  distinct  values 
for  Kp  and  Gp. 

4.4  Numerical  Solutions  of  Equations 

4.4.1  Uniaxial  Stress 

A  PC  based  computer  program  to  solve  the  governing  equations  under  one-dimensional 
stress  conditions  was  developed.  This  program  generates  stress-strain  plots  for  arbitrary  strain 
rate  and  confining  pressure  histories.  A  diagonally  implicit  Runge-Kutta  (DIRK)  numerical 
scheme  described  in  Appendix  D  was  used  for  the  solution  of  three-dimensional  constitutive 
equations  under  uniaxial  stress  conditions.  The  uniaxial  responses  of  AD-85  ceramic  to  different 
strain  rates  and  confining  pressures  were  generated.  The  porous  AD-85  ceramic  was  analyzed 
using  the  proposed  ceramic  model.  The  compressive  stress-strain  curves  generated  by  this  model 
are  shown  in  Figure  4.1.  This  figure  shows  the  effects  of  both  strain  rate  and  confining  pressure. 
The  stress-strain  curves  that  correspond  to  strain  rates  of  100  and  2000/sec.  exhibit  non-linear 
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Figure  4.1.  Stress-Strain  Curves  for  AD-85  Ceramic  at  Different  Strain  Rates  and  Confining 
Pressures. 
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behavior  that  is  due  to  cracking  under  compression.  At  extremely  high  strain  rates,  the  stress 
increases  sufficiently  to  satisfy  the  yield  criterion  and  the  behavior  becomes  dominated  by  plastic 
flow.  The  transition  from  brittle  to  ductile  behavior  as  a  function  of  strain  rate  is  illustrated  in 
Figure  4.2.  The  strength  enhancement  resulting  from  increasing  strain  rates  is  due  to  the  limiting 
crack  velocity  (beginning  at  a  strain  rate  of  about  100/sec.).  If  viscoplastic  flow  is  suppressed, 
the  strength  continues  to  increase  rapidly  with  strain  rate  as  shown  by  the  solid  curve  in  the 
figure.  When  viscoplastic  flow  is  included,  the  peak  stress  cannot  exceed  the  flow  strength  and 
microcracking  becomes  less  important.  Figure  4.3  illustrates  the  transition  from  brittle  to  ductile 
behavior  as  a  function  of  confining  pressure.  Initially,  increased  confining  pressure  delays  the 
onset  of  shear  cracking  and  the  strength  gradually  increases  as  shown  in  the  figure.  This  trend 
continues  until  the  stress  state  reaches  a  threshold  value  for  viscoplastic  pore  collapse.  Once  this 
value  is  attained,  the  strength  becomes  more  dependent  on  the  viscoplastic  properties  than  the 
cracking  properties. 

Figure  4.3  also  includes  data  obtained  from  uniaxial  stress  (confined  and  unconfined)  and 
uniaxial  strain  experiments.  As  shown  in  the  figure,  the  variation  of  strength  as  a  function  of 
pressure  is  successfully  modeled  using  the  microphysical  model.  These  results  are  described  in 
detail  in  Reference  [47]. 

4.4.2  Plate  Impact  Simulation 

In  the  analysis,  plate  impact  experimental  data  were  used  to  calibrate  the  model  constants 
for  AD-85  Alumina  and  for  Titanium  Diboride  [48].  In  these  experiments,  either  the  velocity 
history  was  recorded  using  a  velocity  interferometer  (VISAR)  or  the  stress  history  at  the  interface 
of  the  target  and  back-up  PMMA  was  recorded  using  a  manganin  stress  gauge.  When  a  flyer 
plate  impacts  a  target  plate  at  low  velocities,  fracture  is  induced  in  the  target  by  the  tension 
arising  from  the  interaction  of  reflected  waves  from  the  stress-free  planes.  This  interaction  leads 
to  growth  and  coalescence  of  microcracks,  thus  spalling  the  target.  At  high  impact  speeds,  the 
damage  growth  during  the  initial  compression  of  the  target  alters  the  subsequent  spallation 
process. 

We  considered  three  different  experiments  for  calibrating  and  validating  the  model 
constants.  In  the  first  experiment,  the  target  was  AD-85  Alumina  and  the  impact  velocity  was 
570  m/s,  which  is  above  the  HEL.  Upon  this  high  velocity  impact,  the  shock  stress  in  the  target 
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Figure  4.2.  Effects  of  Dynamic  Cracking  and  Viscoplastic  Flow  on  the  Compressive  Strength 
of  AD-85. 


8 


Shear  Strength  (kbar) 


35 

30 

25 

20 

15 

10 


0 

0  5  10  15  20  25  30  35  40 

Mean  Stress  (kbar) 

Figure  4.3.  Effect  of  Confining  Pressure  on  Strength  of  AD-85. 
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exceeded  the  Hugoniot  elastic  limit  (HEL)  and  the  corresponding  stress  signal  exhibited  a 
ramping  ’plastic’  wave  as  shown  in  Figure  4.4.  Due  to  lack  of  any  physical  evidence,  we  can 
only  speculate  that  either  pore  crushing  or  strain  rate  dependent  plastic  flow  could  generate  such 
a  ramping  wave.  Pore  crushing  will  lead  to  randomly  oriented  microcracks  during  the  initial 
compression.  The  EOS  and  compressive  strength  constants  were  obtained  from  Yaziv  [34]  and 
were  corrected  for  the  10%  voids  in  AD-85  with  the  Mackenzie  formulae.  Figure  4.4  shows  the 
fit  to  the  data  that  was  achieved  with  the  calibrated  model  constants.  Upon  tensile  loading,  the 
microcracks  in  all  orientations  are  sufficiently  large  to  satisfy  the  Griffith  criterion  almost 
simultaneously.  This  leads  to  isotropic  damaged  moduli  behavior  and  rapid  microcrack  growth 
to  a  pulverized  state.  The  model  predicted  an  almost  complete  pulverization  of  the  target  plate 
by  4  ps,  as  verified  by  post  shock  observations. 

In  the  second  experiment,  a  double  flyer  plate  impacted  an  AD-85  Alumina  target  at  a 
lower  velocity  of  293  m/s.  This  impact  generated  shock  stresses  well  below  the  HEL.  The 
VISAR  data,  shown  in  Figure  4.5,  shows  a  spall  signal  after  the  arrival  of  the  release  wave,  and 
then  a  recompaction  wave  signal  from  the  second  flyer.  With  the  isotropic  microcrack  option 
we  could  match  the  spall  signal  but  could  not  fit  the  subsequent  recompaction  signal.  A  few 
simulations  with  this  option  indicated  that  it  was  not  possible  to  match  the  entire  gauge  signal 
over  a  wide  range  of  model  parameters.  The  isotropic  damage  model  overpredicted  the  damage 
level  in  the  target  at  velocities  below  the  HEL.  However,  using  a  model  option  in  which  cracks 
are  allowed  to  coalesce  along  the  perpendicular  to  the  shock  wave  propagation  direction,  we 
achieved  matching  of  the  double  flyer  plate  impact  experimental  data  as  shown  in  Figure  4.5. 

The  third  experiment  was  a  flyer  impact  on  a  TiB2  target  plate  at  a  velocity  of  719  m/s. 
The  manganin  gauge  stress  history  is  compared  with  the  simulation  results  in  Figure  4.6.  The 
experiment  shows  a  descending  compressive  stress  as  a  function  of  time  and  then  a  fairly  level 
spall  signal.  The  growth  of  randomly  oriented  microcracks  in  the  non-porous  TiB2  ceramic 
during  both  compression  and  tension  is  postulated  to  explain  the  trend  in  the  data.  A  fairly  good 
fit  was  obtained  by  using  the  isotropic  microcracking  option.  The  Griffith  criterion  parameters 
were  set  close  to  zero  to  initiate  crack  growth  upon  arrival  of  the  shock  wave. 
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Figure  4.4.  Manganin  Stress  Gauge  Data  and  Model  Simulation  for  Single  Flyer  Impact  on 
Alumina  Plate  at  Velocity  570  m/sec. 
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We  performed  several  plate  impact  tests  on  Silicon  Nitride  (Si3N4)  targets  and 
determined  the  Hugoniot  Elastic  Limit  (HEL)  as  12  GPa.  For  modeling,  we  considered  two 
different  tests  at  impact  velocities  175  m/s  and  800  m/s.  The  stress  amplitude  (3.0  Gpa)  in  the 
low  velocity  test  is  therefore  within  the  elastic  range.  Damage  under  compression  is  not  expected 
to  occur  at  this  stress  level.  The  target  fails  due  to  the  tensile  waves  generated  due  to  wave 
reflections  from  the  stress  free  back  surface  of  the  flyer  and  the  target  plates.  However,  above 
the  HEL  stress  level,  complex  inelastic  deformations  due  to  microplasticity  and  microcracking 
under  compressive  loading  usually  occur  in  the  target  material.  In  fact  the  failure  process  is 
almost  complete  prior  to  the  arrival  of  any  tensile  waves. 

Rajendran  et  al.  [49]  modeled  the  ‘below  HEL’  and  ‘above  HEL’  experiments  using  the 
EPIC-2  code.  Figure  4.7a  compares  the  model  generated  stress  history  with  the  gauge  signal 
from  the  low  velocity  experiment.  The  flat  top  (between  A  and  B )  in  this  figure  indicates  that 
there  is  no  microfracturing  under  compression  and  the  deformation  is  elastic.  However,  the  target 
spalled  due  to  the  tensile  waves  as  indicated  by  the  signal  between  C  and  D.  In  the  absence  of 
spall,  the  stress  history  would  continue  to  decrease,  as  shown  by  the  dotted  line.  In  Figure  4.7b, 
the  model  simulated  stress  history  inside  the  target  is  shown.  The  shock  wave  arrives  at  A  and 
the  material  remains  in  compression  until  the  unloading  wave  arrives  at  C.  Tensile  loading 
begins  at  D,  followed  by  stress  relaxation  due  to  tensile  damage  between  E  and  F.  The  dashed 
line  in  Figure  4.7b  shows  the  damage  evolution  inside  the  target.  Complete  failure  is  assumed 
when  damage  reaches  a  value  of  one. 

In  ceramics,  it  is  easier  to  nucleate  and  propagate  microcracks  under  tension  than  under 
compression.  Therefore,  in  tension,  n+  in  Equation  (74)  is  assumed  to  be  equal  to  one. 
However,  in  compression,  n*  is  usually  small  («  1)  and  is  calibrated  by  matching  the  gauge 
signal. 

Figure  4.8a  compares  the  model  generated  stress  history  with  the  gauge  signal  from  the 
high  velocity  (above  HEL)  experiment.  In  this  figure,  the  stress  drop  between  points  A  and  B 
signifies  compressive  damage  growth  in  the  target.  Figure  4.8b  shows  the  model  generated  stress 
and  damage  histories  inside  the  target.  Compressive  damage  initiation  and  growth  leads  to  stress 
relaxation  between  B  and  C.  At  C,  the  release  wave  from  the  flyer  unloads  the  stress.  In 
general,  the  ability  of  the  model  in  reproducing  the  measured  stress  histories  was  extremely  good. 
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Figure  4.8.  High  Velocity  (800  m/s)  Plate  Impact  Test  with  Si3N4  Target 

(a)  Comparison  Between  Model  and  Gauge  Signal, 

(b)  Model  Generated  Stress  (Solid  Line)  and  Damage  (Dashed  Line)  Histories 
Inside  the  Target. 


4.5  Bar-on-Bar  Impact  Simulation 

The  model  constants  were  calibrated  using  the  data  from  both  plate  impact  and  split 
Hopkinson  bar  experiments.  To  further  evaluate  the  model  constants.  Grove  and  Rajendran  [27] 
simulated  a  rod  impact  experiment,  using  the  EPIC-2  finite  element  code.  The  computed  stress 
history  was  compared  with  the  measured  stress  history  inside  the  target  rod.  The  model  constants 
required  further  adjustment  to  obtain  agreement  with  these  experimental  data. 

An  experiment  was  performed  in  which  a  short  AD-85  ceramic  rod  impacted  a  long  AD- 
85  ceramic  rod  at  a  velocity  of  1 10  m/s.  The  projectile  and  target  rod  lengths  were  50.8  mm  and 

152.4  mm,  respectively,  and  both  rods  were  12.7  mm  in  diameter.  A  manganin  gauge,  embedded 

25.4  mm  from  the  free  end  of  the  target  rod,  was  used  to  measure  the  stress  history  in  the  rod. 
Figure  4.9  illustrates  the  experimental  configuration.  Both  rods  broke  into  small  pieces  as  a 
result  of  the  impact. 

The  rod  impact  experiment  was  simulated  as  an  axisymmetric  problem  using  the  EPIC-2 
finite  element  code.  The  AD-85  ceramic  was  modeled  as  an  elastic-plastic  material  with  10 
percent  porosity.  The  microcrack  model  constants  were  derived  from  split  Ilopkinson  bar  data 
and  a  plate  impact  experiment.  In  the  first  simulation,  damage  due  to  microcracking  was 
suppressed  (the  ceramic  remained  intact  and  elastic)  and  the  results  were  compared  with  the 
experimental  stress  history.  As  can  be  seen  in  Figure  4.10,  the  simulated  stress  history  had  a 
higher  peak  than  the  measured  stress  history.  The  lower  stress  amplitude  in  the  experiment 
indicates  either  degradation  of  the  material  or  some  sort  of  stress  release  due  to  fracture  at  the 
impact  plane. 

In  the  second  simulation,  the  ceramic  was  allowed  to  undergo  damage.  The 
microcracking  was  modeled  using  the  precalibrated  model  constants.  As  Figure  4.11  indicates, 
the  computed  peak  stress  was  significantly  lower  than  the  measured  stress.  The  model  apparently 
overpredicted  the  damage  at  the  impact  end,  and  release  waves  from  the  fractured  end  reduced 
the  amplitude  of  the  initial  shock  wave.  With  lack  of  any  meaningful  data  for  the  microcrack 

propagation  rate,  the  model  parameters  n*,  n^,  nj\  and  n J  (see  Equation  (74))  were  arbitrarily 
assumed  to  be  equal  to  one. 
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Figure  4.9.  Configuration  of  AD-85  Bar-on-Bar  Impact  Experiment. 
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Figure  4.10.  Stress  History  from  Elastic  (No  Microcracking)  Simulation  of  AD-85  Bar-on-Bar 
Impact  Experiment. 
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Figure  4.11.  Model  Generated  Stress  History  with  Microcrack  Growth  Rate  Parameters  (n*, 

n 2 ,  and  n 2 )  Equal  to  1.0  (Excessive  Crack  Growth)  for  AD-85  Bar-on-Bar 

Impact  Experiment. 


It  was  possible  to  match  the  experimental  stress  profile,  as  shown  in  Figure  4.12,  by 
adjusting  the  model  constants  to  limit  the  microcrack  growth  rates  for  both  tensile  (mode  I)  and 
shearing  (modes  II  and  III )  conditions.  The  corresponding  values  for  the  constants  are  given 
in  Table  4.1.  Using  the  new  model  constants,  we  were  still  able  to  match  the  plate  impact  and 
Hopkinson  bar  experimental  results. 


TABLE  4.1 

MODEL  CONSTANTS  FOR  AD-85  CERAMIC 


SYMBOL 

VALUE 

DESCRIPTION 

Kic 

3  MPaJrn 

Static  fracture  toughness 

P 

0.72 

Coefficient  of  friction 

N*o 

1.83  xlO10  m~2 

Microcrack  density  coefficient 

ao 

58  x  10-6  m 

Initial  microcrack  size 

+ 

nl 

1.0 

Microcrack  growth  rate  constants  for  mode  I 

+ 

n2 

0.07 

ni 

0.1 

Microcrack  growth  rate  constants  for  modes  II 
and  III 

n2 

0.07 

4.6  Long  Rod  Penetration  Simulation 

To  verify  the  applicability  of  the  microphysical  based  ceramic  model,  a  hypothetical 
penetration  experiment  was  simulated.  In  the  simulation,  a  long  rod  with  an  aspect  ratio  (L/D, 
L=length,  D=diameter)  of  8  impacts  a  50.8  mm  thick  ceramic  plate  which  is  confined  between 
two  25.4  mm  thick  steel  plates  as  shown  in  Figure  4.13. 

The  EPIC-2  penetration  calculation  was  performed  for  25  microseconds.  During  this 
time,  the  projectile  nearly  penetrated  the  first  steel  plate.  The  main  objective  of  the  analysis  was 
to  understand  the  damage  evolution  inside  the  ceramic  at  various  locations.  Figure  4.14  indicates 
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Figure  4. 12.  Model  Generated  Stress  History  with  New  Set  of  Constants  for  AD-85  Bar-on-Bar 
Impact  Experiment. 
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the  locations  of  the  elements  in  the  ceramic  plate  where  time  history  information  was  obtained 
from  the  simulation.  Element  3665  is  near  the  top  surface  of  the  ceramic  layer  (about  30.5  mm 
from  the  axis  of  symmetry),  element  6402  is  near  the  center  of  the  ceramic  layer  (on  the  axis  of 
symmetry),  and  element  9083  is  near  the  bottom  surface  of  the  ceramic  layer  (about  40.6  mm 
from  the  axis  of  symmetry).  These  three  elements  were  selected  for  the  following  reasons: 

1.  Element  3665  is  expected  to  fracture  completely  due  to  tensile  loading  since  it  is 
away  from  the  axis  of  symmetry. 

2.  Element  6402,  on  the  axis  of  symmetry,  is  expected  to  experience  compressive 
fracture  prior  to  any  tensile  fracture. 

3.  Since  element  9083  is  located  at  the  confined  back  surface  of  the  ceramic  and  away 
from  the  axis,  it  is  expected  to  fracture  due  to  combined  compressive  and  tensile 
loading. 

In  addition  to  time  history  plots,  damage  contour  plots  were  generated  at  8,  10,  12,  and 
20  microseconds.  These  plots,  shown  in  Figure  4.15,  illustrate  the  evolution  of  microcrack 
damage  in  the  ceramic  target  plate.  Damage  initiates  at  locations  near  element  3665  (see  Figure 
4.15a),  well  before  the  penetrator  makes  any  contact  with  the  ceramic  layer.  The  brittle  fracture 
of  the  ceramic  occurs  as  soon  as  it  is  loaded  in  tension  due  to  lateral  release  waves  at  6-8 
microseconds.  The  expanding  region  between  the  damage  contours  in  the  figure  represents  the 
pulverized  ceramic  material.  The  damage  is  progressive  as  seen  by  the  contours  movement  in 
Figure  4.15.  A  small  portion  (region  A  in  Figure  4.15c)  of  the  ceramic  directly  under  the 
penetrator  remains  intact  while  the  damage  that  initiated  near  element  3665  propagates  well  into 
the  ceramic.  At  20  microseconds,  the  strong  compressive  loading,  induced  by  the  approaching 
rod,  completes  the  fracturing  process.  As  Figure  4.15d  indicates,  the  ceramic  material  directly 
beneath  the  penetrating  rod  pulverizes  completely  before  the  rod  reaches  the  ceramic.  A  mesh 
plot  of  the  fractured  ceramic  is  shown  in  Figure  4.16;  the  dark  region  is  the  pulverized  ceramic. 

The  time  history  plots  in  Figures  4.17  through  4.19  provide  information  concerning 
elements  3665,  6402,  and  9083.  These  figures  include  maximum  principal  stress,  pressure,  and 
damage  time-history  plots.  Elements  3665  and  6402  experience  only  compressive  pressures  as 
can  be  seen  from  Figure  4.17a  and  4.18a.  The  compressive  pressure  is  positive  in  the  plots. 
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Microcrack  Damage  Contours 
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Microcrack  Damage  Contours  at  (a)  8  ps,  (b)  10  ps,  (c)  12  ps,  and  (d)  20  ps 
(concluded). 
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Figure  4.17.  Time  History  Plots  for  Element  3665:  (a)  Pressure,  and  (b)  Maximum  Principal 
Stress  (Solid  Line)  and  Damage  (Dashed  Line). 
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Stress  (Solid  Line)  and  Damage  (Dashed  Line). 


However,  element  9083  experiences  a  small  tensile  pressure  at  approximately  t=13  microseconds, 
as  shown  in  Figure  4.19. 

Plots  of  principal  stress  and  damage  versus  time  are  presented  in  Figures  4.17b,  4.18b, 
and  4.19b  for  elements  3665,  6402,  and  9083,  respectively.  The  model  restricts  the  maximum 
positive  (tensile)  principal  stress  to  a  few  kilobars  (<3  kbars).  This  is  achieved  by  the  model 
through  stress  relaxation  due  to  damage  evolution.  Element  3665  experienced  only  Mode  I 
microcrack  growth,  while  the  microcracking  in  element  6402  began  in  Modes  II  and  III  finished 
in  Mode  I.  However,  the  crack  growth  under  compression  (shear  cracking)  is  relatively  small 
compared  to  the  tensile  growth  and  is  not  discernable  in  the  figure.  Element  9083,  like  3665, 
experienced  only  rapid  Mode  I  microcrack  growth.  These  time  history  plots  indicate  that  the 
model  predicted  damage  evolution  in  the  ceramic  according  to  the  generalized  Griffith  criterion 
and  stress  relaxation  according  to  micromechanics.  The  stress  relaxation  under  tensile  loading 
is  clearly  demonstrated  in  the  simulation.  Unlike  models  in  which  damage  does  not  influence 
the  strength,  this  advanced  ceramic  model  incorporates  the  effects  of  damage  on  both  strength 
and  stiffness.  Thus,  the  use  of  this  newly  developed  ceramic  model  in  penetration-into-ceramic 
calculations  should  provide  hydrocodes  with  the  necessary  predictive  capabilities. 
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Section  5 

Summary  and  Recommendations 


The  development  of  advanced  constitutive  models  for  use  in  hydrocode  calculations  was 
reported.  A  ductile  failure  model  (the  RDG  model)  for  metals  and  a  microcracking  model  for 
ceramics  were  presented.  In  the  model  formulations,  the  evolution  equations  were  directly  related 
to  the  physics  operating  at  the  microscopic  level.  One  goal  in  the  model  development  was  to 
describe  the  physical  process  with  meaningful  mathematical  expressions.  Another  important 
philosophy  was  to  keep  the  model  constants  to  a  minimum. 

5.1  RDG  Model  Summary 

The  RDG  model  is  a  continuum  mechanics  based  ductile  failure  model.  The  model 
formulation  is  three-dimensional  and  therefore  applicable  under  a  general  stress-strain  state.  The 
RDG  model  describes  several  fundamental  aspects  of  the  failure  process.  When  the  failure 
process  is  initiated  due  to  void  nucleation  in  the  material,  the  plastic  flow  becomes  pressure 
dependent.  Also,  under  high  velocity  impact  loading,  most  metals  exhibit  strain  rate  dependent 
behavior.  The  RDG  model  includes  the  effects  of  pressure  and  strain  rate  in  the  constitutive 
(strength)  model  formulation.  An  important  feature  of  the  RDG  model  is  its  accurate  modeling 
of  the  intact  material  using  appropriate  viscoplastic  equations.  The  Bodner-Partom  (BP) 
equations  were  selected  for  this  purpose  although  we  have  shown  that  the  Johnson-Cook  (JC) 
equations  may  also  be  used.  The  state  variables  in  the  RDG  model  formulation  are  the 
equivalent  plastic  strain  rate  in  the  matrix  material  and  the  BP  model  state  variable  Z,  which 
describes  the  loading  history  effects. 

A  pressure  dependent  yield  function  serves  as  the  plastic  potential  in  the  derivation  of 
plastic  strain  rate  for  the  void  contained  aggregate.  Since  the  aggregate  plastic  strain  rates  are 
directly  related  to  the  equivalent  plastic  strain  rate  of  the  matrix  material,  the  strain  rate  and 
loading  history  effects  enter  into  the  failure  model  formulation  through  the  Bodner-Partom  model. 
The  RDG  failure  model  does  not  require  separate  loading  or  unloading  conditions.  Both  elastic 
and  plastic  rate  components  are  taken  to  be  always  nonzero,  and  the  same  relations  are  intended 
to  hold  under  loading  and  unloading  conditions  as  in  the  Bodner-Partom  model.  When  the 
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equivalent  plastic  strain  rates  are  numerically  very  small  (<0.000 1/sec.),  the  deformation  is 
considered  to  be  elastic. 

Another  appealing  feature  of  the  RDG  failure  model  is  the  void  growth  law.  Assuming 
that  the  void  volume  fraction  growth  rate  is  equal  to  the  summation  of  the  normal  plastic  strain 
rate  components  (dilatation)  removes  the  requirement  for  a  new  evolution  law  for  the  void  growth 
rate.  This  modeling  approach  significantly  reduces  the  number  of  model  constants.  While  the 
Bodner-Partom  model  usually  requires  five  constants,  the  determination  of  those  constants  is 
fairly  straight  forward  and  Rajendran  et  al.  [32]  reported  a  procedure  to  calibrate  the  constants 
using  data  from  a  limited  number  of  standard  impact  tests. 

In  the  RDG  model,  there  are  three  constants  for  void  nucleation,  and  two  constants  that 
control  the  growth  of  voids.  Model  constants  can  be  determined  from  numerical  simulations  of 
a  plate  impact  experiment.  The  initial  value  for  the  void  volume  fraction  and  the  threshold 
nucleation  stress  can  be  adjusted  to  match  the  time  of  arrival  of  the  experimental  spall  signal. 
The  standard  deviation  of  the  nucleation  stress  distribution  is  arbitrarily  set  to  one  third  or  one 
fourth  of  this  threshold  stress.  The  initial  slope  and  amplitude  of  the  spall  signal  guides  the 
selection  of  the  void  growth  (or  yield  function)  constants. 

The  RDG  model  was  incorporated  into  the  EPIC-2  finite  element  code.  A  numerically 
stable  and  accurate  solution  technique  based  on  a  diagonally  implicit  Runge-Kutta  (DIRK) 
method  was  used  to  solve  the  complex  and  stiff  governing  differential  equations.  Various  aspects 
of  the  numerical  scheme  were  demonstrated  in  References  [1]  and  [50].  For  solving  only  the 
Bodner-Partom  equations,  an  iterative  radial  return  scheme  has  been  successfully  implemented. 
The  BP  model  constants  had  already  been  determined  for  these  materials  and  reported  by 
Rajendran  et  al.  [51].  The  RDG  model  accurately  matched  several  of  the  plate  impact 
experiments.  RDG  model  constants  for  OFHC  copper,  tantalum,  Armco  iron,  and  1020,  MAR 
200,  MAR  250,  AF1410,  C1008,  and  HY100  steels  were  determined  and  are  presented  in  this 
report. 

The  capability  of  the  RDG  model  was  extended  to  multiple  impact  loading  conditions. 
The  model  accurately  predicted  the  failure  process  in  a  twice-impacted  target  plate.  Both  void 
growth  and  collapse  were  modeled  in  this  case.  To  describe  the  pore  collapse  correctly,  a  critical 
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value  of  void  volume  fraction  was  calibrated  for  collapse.  For  complete  spallation,  the  critical 
void  volume  fraction  was  arbitrarily  assumed  to  be  1.0  for  most  metals. 

We  also  successfully  modeled  spallation  under  a  three-dimensional  strain  state  using  the 
RDG  model.  The  model  constants  were  determined  from  a  one-dimensional  (strain)  plane  plate 
impact  test.  Using  those  constants,  the  RDG  model  reproduced  multiple  spall  regions  in  a  cone 
target  under  a  three-dimensional  stress-strain  state.  This  demonstration  greatly  increased  our 
confidence  in  the  generality  of  this  model,  even  when  calibrated  with  one-dimensional 
experiments.  We  also  established  that  the  correct  spall  patterns  in  the  cone  can  only  be  obtained 
by  using  a  void  nucleation  and  growth  model  such  as  the  RDG  model. 

5.2  Ceramic  Model  Summary 

Initially,  a  relatively  simple  fragmentation  based  ceramic  model  (the  modified  TCK 
model)  was  developed  and  reported  in  Reference  [11].  However,  based  on  the  complexity  of 
impact  damage  evolution  under  both  compressive  and  tensile  loading,  an  advanced  microphysical 
model  was  developed.  This  model  is  capable  of  modeling  inelastic  strains  due  to  pore  collapse, 
plastic  flow,  and  microcracking  in  ceramics.  This  advanced  model  has  been  incorporated  into 
the  EPIC-2  computer  code.  Plate  impact,  rod-on-rod  impact,  and  long  rod  penetration  into 
ceramic  target  experiments  were  simulated.  In  the  simulation,  when  the  ceramic  material 
completely  pulverizes  at  impact  velocities  above  the  HEL,  the  isotropic  damage  model  is 
sufficient  to  describe  the  tensile  failure.  The  pore  crushing  and  strain  rate  dependent  plasticity 
under  compression  are  modeled  using  the  state  variable  and  pressure  dependent  plastic  flow 
theories. 

The  EPIC-2  code  simulation  of  a  plate  impact  test  on  AD85  ceramic  using  the  ceramic 
model  reproduced  the  measured  stress-time  history.  Since  this  is  a  porous  ceramic,  the  inelastic 
("plastic")  ramping  of  the  stress  was  modeled  through  the  plastic  pore  collapse  option  of  the 
ceramic  model.  The  damage  in  AD-85  and  TiB2  at  impact  velocities  above  HEL  was  modeled 
by  treating  the  damage  as  a  randomly  oriented  scaler  variable. 

Since  microdamage  occurs  at  shock  levels  well  below  the  HEL,  it  is  important  to  model 
fracture  initiation  under  compression.  Using  the  microphysical  model  which  contains  this  feature, 
a  plate  impact  test  on  a  silicon  nitride  target,  at  a  velocity  well  below  the  HEL,  was  successfully 
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modeled.  The  model  reproduced  the  spall  signals  extremely  well.  We  successfully  modeled 
impact  damage  evolution  in  silicon  nitride  under  both  tensile  and  compressive  loading  conditions. 
The  strength  and  stiffness  degradations  were  modeled  through  an  internal  state  variable  based 
scaler  damage  parameter  y  representing  the  crack  density  in  the  material.  We  demonstrated  the 
capabilities  of  the  microphysical  model  by  modeling  both  ‘below  HEL’  and  ‘above  HEL’  plate 
impact  experiments. 

A  long  rod  penetration  into  a  confined  ceramic  target  was  also  simulated  using  the  EPIC- 
2  code.  The  impact  behavior  of  the  ceramic  target  was  described  by  the  advanced  ceramic 
model.  The  code  analysis  showed  the  details  of  damage  evolution  under  multi  axial  loading. 
Physically  meaningful  results  were  obtained  from  the  code  simulation.  The  capability  of 
obtaining  realistic  solutions  from  the  hydrocode  calculations  using  such  physically  based  models 
is  an  important  advantage  in  advanced  impact  design  feasibility  studies. 

5.3  Recommendations 

5.3.1  Metals 

The  results  from  the  analysis  of  the  application  problems  suggest  the  need  for  additional 
in-depth  analyses  supported  by  experimental  data.  For  example,  to  substantiate  the  tensile 
necking  analysis,  high  speed  photographic  split  Hopkinson  bar  tests  on  shallow  notch  specimens 
should  be  conducted.  The  experimentally  determined  necking  profiles  should  be  compared  with 
the  RDG  model  predictions.  For  this  purpose,  selection  of  a  material  in  which  microvoids  have 
been  observed  under  dynamic  tensile  loading  conditions  is  recommended. 

It  is  desirable  to  perform  instrumented  penetration  experiments  for  validating  the  RDG 
model  prediction.  Radiographs  of  the  projectile  and  target  during  the  penetration,  and  post-failure 
analysis  on  the  target  would  identify  additional  model  requirements.  These  radiographs  are 
particularly  important  for  verifying  the  spall  fragmentation  studies. 

5.3.2  Ceramics 

The  generality  of  the  new  ceramic  model  requires  additional  validation.  Simulations  of 
plate  impact  tests  with  different  options,  such  as  (1)  elastic-cracking,  (2)  elastic-viscoplastic 
without  pore  collapse,  (3)  elastic-viscoplastic  with  pore  collapse,  and  (4)  elastic-cracking- 
viscoplastic,  should  be  performed.  Since  the  impact  response  of  ceramics  differs  significantly 
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from  ceramic  to  ceramic,  it  is  important  that  other  ceramic  materials,  such  as  boron  carbide, 
silicon  carbide,  aluminum  nitride,  etc.,  are  also  considered  in  the  model  evaluation.  The  effects 
of  pulverization  on  the  stress  wave  profiles  must  be  investigated. 

Additional  efforts  are  required  for  the  bar-on-bar  and  long  rod  penetration  experiments. 
Here  again,  the  various  options  in  the  ceramic  model  should  be  invoked  and  the  corresponding 
simulation  results  should  be  compared.  A  model  constants  evaluation  scheme  needs  to  be  firmly 
established.  The  development  of  a  standard  procedure  is  essential  so  that  material  constants  for 
various  ceramics  can  be  systematically  derived. 

The  ability  of  the  ceramic  model  to  predict  the  penetration  depth  of  a  long  rod 
penetrating  both  confined  and  unconfined  ceramics  must  be  established.  This  will  require  data 
from  diagnostic  penetration  experiments.  In  summary,  further  research  is  required  to  validate  the 
newly  developed  microphysical  model. 
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Appendix  A 

Split  Hopkinson  Bar  and  Plate  Impact  Data 


This  section  comprises  a  compendium  of  experimental  data  obtained  during  the  course 
of  this  investigation.  The  SHB  data  are  presented  in  terms  of  plots  of  stress  versus  log  strain  rate 
for  a  given  value  of  strain  for  several  metals.  The  plate  impact  data  are  in  terms  of  plots  of 
measured  stress  versus  time  for  several  metals.  Brief  descriptions  of  both  these  test  techniques 
are  also  presented  in  this  section. 

A.1  Split  Hopkinson  Bar  (SHB) 

The  SHB  apparatus  consists  of  a  striker  bar  that  is  made  to  impact  end-to-end  onto  the 
pressure  bar.  The  impact  produces  a  stress  wave  in  the  pressure  bar  whose  duration  is  twice  the 
acoustic  transit  time  of  the  striker  bar  (0.3  ms  in  the  bar  used  at  the  University  of  Dayton 
Research  Institute).  The  specimen  is  placed  at  the  other  end  of  the  pressure  bar  (bar  1)  as  shown 
in  Figure  A. la.  The  specimen  is  also  connected  to  the  transmission  bar  (bar  2).  The  striker, 
pressure,  and  transmission  bars  are  all  the  same  material  and  diameter  (12.7-mm  diameter  Inconel 
in  these  experiments).  The  stress  wave  that  travels  down  the  pressure  bar  is  partially  transmitted 
and  partially  reflected  by  the  specimen.  The  reflected  and  transmitted  waves  in  the  bars  are 
detected  by  strain  gauges.  Analysis  of  the  strain-gauge  signals  yields  the  load  and  displacement 
history  of  the  specimen.  Detailed  discussions  on  the  strain-gauge  analysis  can  be  found  in 
References  [A.l]  and  [A.2]. 

For  tensile  measurements,  a  collar  is  placed  around  the  specimen,  and  the  specimen  is 
screwed  into  the  pressure  and  transmission  bars.  (See  Figure  A. lb.)  The  initial  compression 
wave  is  transmitted  through  the  collar,  and  reflects  as  a  tensile  wave.  When  the  tensile  wave 
returns  to  the  specimen,  the  collar  falls  away  and  the  specimen  sustains  the  tensile  load.  The 
data  from  the  Hopkinson-bar  strain  gauges  are  used  to  calculate  engineering  stress  and  strain  in 
the  sample,  which  are  converted  to  true  stress  and  strain  until  the  onset  of  localized  deformation 
(necking)  by  using  the  conventional  relationships. 
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A.2  Plate  Impact  Experiments 


Many  discussions  of  planar  impact  loading  are  available  (References  [A.3]  and  [A.4]). 
Plate  impact  tests  provide  a  loading  path  that  is  very  different  from  conventional  SHB 
compression  or  tension  tests.  The  deformation  is  that  of  one-dimensional  strain,  and  the  mean 
stress  is  generally  very  high  compared  to  that  in  SHB  tests.  Strain  rates  are  105  s'1  or  higher. 
The  material  undergoes  compression  immediately  followed  by  tension.  Thus,  plate  impact 
experiments  are  essential  for  calibrating  and  validating  high  strain  rate  material  models  that  aspire 
to  general  applicability.  Specifically,  plate  impact  data  may  be  interpreted  to  infer  compression 
and  tensile  yield  strengths  and  failure  parameters. 

The  plate  impact  tests  were  conducted  with  three  objectives:  (1)  determination  of 
Hugoniot  Elastic  Limit,  (2)  determination  of  the  unloading  path  from  the  free  surface  velocity 
history,  and  (3)  determination  of  the  threshold  conditions  associated  with  onset  of  spall  fracture. 

Impact  induces  an  elastic  shock  and  a  plastic  shock  in  the  target.  The  amplitude  of  the 
elastic  shock  is  The  Hugoniot  elastic  limit,  G^pj ,  is  the  maximum  stress  for  one¬ 

dimensional  elastic  wave  propagation.  This  stress  level  is  a  material  property,  and  above  this 
level  the  material  flows  plastically.  The  stress  (J^  can  be  determined  from  the  experimentally 
obtained  free  surface  velocity  of  the  target  that  corresponds  to  the  elastic  shock,  u^j . 

aHEL  =  PCL  UHEL  (A.l) 


where  cL  is  the  elastic  sound  speed.  However,  when  a  stress  gauge  is  used  to  measure  GhpT  we 
need  the  impedance  match  solution  to  interpret  the  data.  Since  the  measured  stress  is  the  stress 
that  is  transmitted  into  the  PMMA  wedge  window  (see  Figure  A.2),  this  G^gL  of  the  PMMA 
must  be  multiplied  by  a  factor  in  order  to  obtain  the  ghfi  of  the  target  material.  The 
corresponding  impedance  relationship  is. 


CHEL  ~ 


7  +  7  '\ 

j  rp  ^  £j  p  * 


Tz, 


>HEL 


(A.2) 


where  Zj  and  Zp  are  the  impedance  (pc)  of  the  target  and  PMMA,  respectively.  In  our 
experiments,  Zp  =  33.5  X  105  Pa-sec/m. 
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A.3  Materials  Data 


The  plots  of  strength  variation  with  respect  to  strain  rate  for  several  metals  are  given  in 
Figures  A.3  -  A.25.  This  data  is  presented  in  terms  of  effective  stress  versus  strain  rate  for  a 
given  strain  level.  This  data  is  generated  from  quasi-static  and  Split  Hopkinson  bar  terms. 

The  measured  stress  versus  time  histories  from  plate  impact  tests  for  several  metals  and 
ceramics  are  plotted  in  Figures  A. 26  -  A.45.  This  data  is  valuable  for  material  model  validation. 
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Figure  A-3  Stress  (strength)  versus  strain  rate  f  >r  MAR-M  300  Steel 
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Figure  A-4  Stress  (strength)  versus  strain  rate  for  3003  Aluminum 
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versus  strain  rate  for  Inconel  600 


Stress  (strength)  versus  strain  rate  for  Inconel  718 
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Figure  A-7  Stress  (strength)  versus  strain  rate  for  AF1410  Steel 
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Stress  (strength)  versus  strain  rate  for  Vanadium 
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Stress  (strength)  versus  strain  rate  for  7017  Aluminum 
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Figure  A- 12  Stress  (strength)  versus  strain  rate  for  1044  Steel 
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(strength)  versus  strain  rate  for  7039-T64  Aluminum 
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Figure  A- 14  Stress  (strength)  versus  strain  rate  for  OFHC  Copper 
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Figure  A- 15  Stress  (strength)  versus  strain  rate  for  HY100  Steel 
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Figure  A- 18  Stress  (strength)  versus  strain  rate  for  Plate  A  and  B,  DU  1  and  DU2 
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Figure  A-20  Stress  (strength)  versus  strain  rate  for  Tungsten 
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Figure  A-21  Stress  (strength)  versus  strain  rate  for  Nickel 
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Figure  A-22  Stress  (strength)  versus  strain  rate  for  Inconel  625  and  pure  Tantalum 
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Figure  A-23  Stress  (strength)  versus  strain  rate  for  52125  Steel  and  Armco  Iron 
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Figure  A-24  Stress  (strength)  versus  strain  rate  for  D-6AC  Steel 


Figure  A-25  Stress  (strength)  versus  strain  rate  for  MAR-M  200  Steel 
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Figure  A-26  Gage  stress-time  history  for  Tungsten  (90%)  plate  impact  test 
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Figure  A-27  Gage  stress-time  history  for  1215  Steel  plate  impact  test 
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Figure  A-28  Gage  stress-time  history  for  4340  Steel  plate  impact  test 
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Figure  A-29  Gage  stress-time  history  for  46100  Steel  plate  impact  test 
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Figure  A-30  Gage  stress-time  history  for  MAR  200  Steel  plate  impact  test 
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Figure  A-31  Gage  stress-time  history  for  MAR  250  Steel  plate  impact  test 


A-33 


STRESS  ( k  bon) 


T I ME  (  microsec) 


Figure  A-32  Gage  stress-time  history  for  Tungsten  (98%)  plate  impact  test 
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Figure  A-33  Gage  stress-time  history  for  Tantalum  plate  impact  test 
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Figure  A-35  Gage  stress-time  history  for  1020  Steel  plate  impact  test 
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Figure  A-36  Gage  stress-time  history  for  Armco  Iron  plate  impact  test 
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Figure  A-37  Gage  stress-time  history  for  C1008  Steel  plate  impact  test 
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Figure  A-39  Gage  stress-time  history  for  AD-90  plate  impact  test 
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Figure  A-40  Gage  stress-time  history  for  AD-995  plate  impact  test 
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Figure  A-41  Gage  stress-time  history  for  Boron  carbide  plate  impact  test 
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Figure  A -42  Gage  stress-dme  h.s.ory  for  Titanium  diboride  p,aK  jmpac|  ^ 
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Figure  A-43  Gage  stress-time  history  for  Silicon  nitride  plate  impact  test 
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Figure  A-44  Gage  stress-time  history  for  Silicon  nitride  plate  impact  test 
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Figure  A-45  Gage  stress-time  history  for  Silicon  nitride  plate  impact  test 
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Appendix  B 

Johnson-Cook  and  Bodner-Partom  Models 


B.l  Johnson-Cook  Model 

The  Johnson-Cook  constitutive  model  [B.l]  provides  realistic  solutions  to  hydrocode 
simulations  of  a  broad  class  of  applications  to  dynamic  events  such  as  impact,  penetration,  and 
explosive  acceleration  of  metals.  The  Johnson-Cook  model  has  the  form: 

y  =  [a  +  BFn]  [l  +  C  Jnfi*]  [l  -  T*m]  (B.l) 


T* 


T  —  Tr0om 
^mlc  ~  ^room 


(B.2) 


where:  Y  is  the  flow  strength,  e  is  effective  plastic  strain,  C*  is  non-dimensional  (effective 
plastic)  strain  rate  (normalized  by  1/sec),  T*  is  homologous  temperature,  Tmll  is  the  temperature 
at  melting,  and  T  is  the  applied  temperature.  The  five  material  constants  are  defined  as  follows: 
A  is  the  static  yield  strength,  B  is  the  work  hardening  coefficient,  n  is  the  work  hardening 
exponent,  C  is  the  strain  rate  coefficient,  and  M  is  the  thermal  softening  exponent.  In  Table  B.l, 
the  values  of  these  constants  are  tabulated  for  various  metals. 

The  JC  model  solution  scheme  in  the  EPIC  code  is  based  on  a  radial  return  method.  The 
model  constants  A,  B,  and  n  are  determined  from  quasi-static  stress-strain  data  either  from  tensile 
or  compressive  tests.  The  strain  rate  dependent  constant  C  is  determined  from  the  slope  of  the 
stress  vs.  strain  rate  plot.  The  temperature  constant  M  is  estimated  from  the  stress  vs. 
temperature  plot.  In  the  code  analysis,  for  a  given  strain,  strain  rate,  and  temperature,  the  von 

Mises  yield  radius  is  calculated  from  the  JC  equation.  The  von  Mises  stress,  ^3  J2  is  given  by, 


\J^J2 


Sij 


'ij 


(B.3) 


where  the  J2  is  the  second  invariant  of  the  stress  deviators,  Sy. 
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not  provided 


B.2  Bodner-Partom  Model 


We  present  the  major  features  of  the  Bodner-Partom  model  [B.2]  here,  as  well  as  an 
extension  of  the  model  for  extreme  work  hardening  materials,  as  proposed  earlier  by  Bodner  and 
Merzer  [B.3]. 

The  total  strain  rate  is  assumed  to  be  decomposable  into  elastic  and  inelastic  components. 


'13 


=  t  -jj  +  t 


1J 


(B.4) 


Both  are  nonzero  for  all  loading/unloading  conditions.  The  elastic  strain  rate  is  related  to  the 
stress  rate  by  the  elastic  constants  (Hooke’s  Law).  The  inelastic  strain  rate  is  assumed  to  be  a 
function  of  stress,  a^,  a  state  variable,  Z,  and  to  follow  the  Prandtl-Reuss  flow  rule 

tPij  =  efj  =  X5ij-  (B.5) 


where  efj  are  the  deviatoric  plastic  strain  rates.  Squaring  Equation  (B.5)  gives 


Df  =  X2J2 


(B.6) 


where  Df  is  the  second  invariant  of  the  plastic  strain  rate. 


D?  =  D2d  exp 


\n 


n  +  1 
n 


(B.7) 


where  D0  is  the  limiting  value  of  the  plastic  strain  rate  in  shear,  and  n  is  a  parameter  that  is 
mainly  related  to  strain  rate  sensitivity.  Z  is  a  measure  of  the  overall  resistance  of  the  material 
to  plastic  flow,  and  it  depends  on  the  loading  history. 


The  plastic  strain  rate  can  be  expressed  in  the  following  form  by  regrouping  Equations 
(B.5)  through  (B.7): 


(B.  8) 
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It  is  assumed  that  Z  is  a  function  of  the  plastic  work,  Wp: 

Z  =  m(z1  -  Z)  Wp 


where  Wp=aiji^j.  Integration  of  Equation  (B.9)  leads  to: 


Z  =  zx 


-(Zi-Zo) 


e 


-mWp 


(B.10) 


Z0  is  the  initial  value  of  Z,  Zx  is  the  maximum  value  that  Z  can  attain,  and  m  is  a  parameter  that 
describes  the  strain  hardening  behavior  of  the  material. 

To  describe  extreme  strain  hardening  material  behavior,  m  can  be  assumed  to  take  the 
following  form: 


m 


ma  +  m1e 


-a 


(B.ll) 


The  expression  for  Z  is  obtained  by  combining  Equations  (B.9)  and  (B.ll)  and  integrating. 


( w0*m1  -w) 

- - - 


+  m0wp 


(B.12) 


The  additional  constants,  m0  and  rrij,  replace  the  original  m.  We  also  need  a  to  model  the 
extreme  strain  hardening  behavior  of  a  material  like  annealed  OFHC  copper  [B.4]. 

The  effect  of  temperature  on  the  flow  stress  is  included  through  the  parameter  n.  Based 
on  quasi-static  experimental  evidence,  the  parameter  n  is  assumed  to  vary  as  an  inverse  function 
of  absolute  temperature  (T): 

n  =  A  +  —  (B.13) 

T 

where  A  and  B  are  model  parameters  and  T  is  expressed  in  either  degrees  Rankine  or  degrees 
Kelvin.  Thus  the  BP  model  contains  five  principal  material  parameters  D0,  Z0,  Zj,  m,  and  n  that 
have  to  be  evaluated  from  high  strain-rate  experiments  under  room  temperature  for  most  metals. 
The  two  secondary  work  hardening  constants  a  and  irq  are  used  to  describe  highly  strain 
hardening  materials  like  OFHC  copper.  In  high  strain  rate  applications,  the  limiting  strain  rate 
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D0  has  been  assumed  to  be  108/sec  for  metals.  This  eliminates  the  task  of  determining  one  more 
model  constant.  Therefore,  the  BP  model  generally  requires  only  4  material  model  constants  for 
room  temperature  applications.  However,  for  high  temperature  applications,  one  additional 
constant  is  required. 

Rajendran  and  Geers  [B.5]  developed  a  user  interactive  program  called  BPSOLVE  for 
the  determination  of  the  first  four  constants  described  above,  based  on  at  least  three  tensile  tests 
at  different  strain  rates.  The  temperature  related  material  constants  can  be  determined  iteratively 
using  high  temperature  SHB  test  data.  Rajendran  et  al.  [B.6,B.7]  determined  model  parameters 
for  several  metals  using  the  BPSOLVE  program.  The  corresponding  model  constants  are  given 
in  Table  B.2. 

In  the  BP  model,  the  plastic  strain  rates  are  expressed  as  a  function  of  stress  components, 
state  variable  Z,  and  the  second  invariant  of  the  stress  deviator.  Therefore,  expressing  the 
strength  Y  in  terms  of  strain  rate  and  state  variables  (as  in  the  JC  model)  is  not  a  straightforward 
task.  In  hydrocodes,  calculation  of  the  deviatoric  stresses  requires  an  iterative  scheme  when 
implementing  the  BP  model  (or  similar  constitutive  relationships).  Initially,  the  plastic  strain 
rates  are  estimated  using  the  deviatoric  stresses  calculated  from  elastic  equations  as, 

-  2G  (e2J  -  .if,)  (B'14) 

The  plastic  strain  rate  expressions  [B.8]  of  the  Bodner-Partom  model  make  these  ordinary 
differential  equations  stiff  under  certain  conditions.  Computation  of  the  deviatoric  stresses 
requires  plastic  strain  rate  estimates,  while  computation  of  the  plastic  strain  rates  requires 
deviatoric  stress  estimates.  Therefore,  the  stress  calculation  requires  an  iterative  scheme  in  the 
numerical  implementation  of  the  BP  model.  Several  articles  on  the  results  from  using  the  BP 
model  in  hydrocodes  can  be  found  in  References  [B.8]  through  [B.13]. 
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TABLE  B.2 

BODNER-PARTOM  MODEL  CONSTANTS 


Material 


C1008 

Steel 


HY100 

Steel 


1020 

Steel 


MAR-200 

Steel 


Armco 

Iron 


OFHC 

Copper 


6061-T6 

Aluminum 


7039-T64 

Aluminum 


Pure 

Tantalum 


W-2 

Tungsten 


Nickel 

200 


MAR-250 

Steel 


AF1410 

Steel 


Zo 

(GPa) 


a 

GPa1 


0.245  46 


NA  NA 


NA  NA 


NA  NA 


NA  NA 


1500  NA  NA 


-2.86  2343 


NA  NA 


NA  NA 


0  0.166  134 


NA  NA 


NA  NA 


NA  NA 


NA  --  The  high  temperature  constants  are  "Not  Available" 
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Appendix  C 

Solution  Scheme 


C.l  Numerical  Solution  of  the  Governing  Equations  in  Sections  3  and  4 

The  governing  equations  described  in  Sections  3  and  4  are  rearranged  and  numerically 
integrated.  A  detailed  description  can  be  found  in  Reference  [C.l].  For  stability,  accuracy,  and 
uniqueness  of  the  solution,  a  second  order  diagonally  implicit  Runge-Kutta  (DIRK)  method  was 
considered  for  integrating  the  system  of  stiff  ordinary  differential  equations.  The  DIRK  scheme 
is  outlined  in  Appendix  D.  The  DIRK  variables  are  integrated  using  the  estimated  rates  for  the 
current  time. 

The  EPIC-2  code  calculates  the  current  strain  rates  in  each  material  element.  Before  any 
void  nucleation  occurs,  the  volumetric  strains  of  the  matrix  and  the  aggregate  are  the  same. 
However,  once  void  nucleation  occurs  in  an  element,  the  incremental  volumetric  strain  of  the 

aggregate,  t  v,  is  the  sum  of  the  incremental  elastic  volume  change,  t  ®,  in  the  aggregate  (in  the 
same  sense  as  Mackenzie’s)  and  the  incremental  volume  change  due  to  the  viscoplastic  growth 
of  voids,  tpii .  This  relationship  can  be  related  directly  to  the  aggregate  elastic  compressible 
strain,  pag,  from  Equation  (13)  as  follows.  The  relationship  between  the  true  volumetric  strain 
and  the  compressible  strain  is,  by  definition, 

1  +|j,  =  e~tv  (C.l) 


The  relationship  between  the  true  plastic  dilatation  and  the  void  volume  fraction  is  given  by  the 
direct  integration  of  Equation  (9)  as. 


1-f  = 


(C.2) 


where  the  plastic  dilatation,  e^,  is  zero  at  f  =  0.  Substituting  Equations  (C.l)  and  (C.2)  into 
Equation  (13),  the  following  relationships  were  obtained: 


\iag  =  e-ev-l 


(C.3) 


C-l 


and 


(C.4) 


e  p 

eV~  ev~e  ii  • 


where  e  ®  is  the  true  elastic  dilatation  strain  of  the  aggregate.  This  development  suggests  that 

the  plastic  strain  rates,  e *?  •,  should  be  numerically  integrated  instead  of  integrating  Equation  (4). 

The  aggregate  pressure,  P,  is  evaluated  by  using  Equation  (C.3)  and  the  procedure  outlined  in 
Section  3.  The  true  elastic  deviatoric  strains  of  the  aggregate  are  given  by  the  equation, 


- 

(  \ 

“ 

.P 

ekk 

- 

- 

eij  ~  eij 


The  aggregate  deviatoric  stresses  can  then  be  determined  from: 

Sij  =  2Geeij  , 

where  G  is  the  degraded  shear  modulus,  as  defined  by  Equation  (11). 


(C .  5) 


(C .  6 ) 


One  numerical  bonus  of  the  above  approach  is  that  the  use  of  Equation  (C.2)  to  calculate 
f  avoids  numerical  integration  of  the  stiff  differential  Equation  (9).  Equations  (4),  (B.8),  and 
(B.9)  lead  to  stiff  differential  equations;  therefore,  the  DIRK  scheme  permits  extremely  small 
time  steps  to  maintain  solution  stability  and  accuracy.  For  numerical  efficiency  (i.e.  to  avoid 
extremely  small  time  steps).  Equations  (B.9)  and  (6)  were  analytically  integrated  to  derive  the 
following  equations  for  the  BP  model  loading  history  parameter  (Z)  and  the  void  volume  fraction 
due  to  nucleation  of  voids  (fn): 


2  ~  Z1  (Z1  Zo)e 


(m0  +  m1  -m) 


a 


(C .  7 ) 


fnm 


(V 


erf 


P+Ym-GN 


erf 


s1 

J 

\ 

~  eN 

-  erf 

'  -e„ 

fi  s2  j 

fi  s; 

(C.8) 


Using  Equation  (C.8)  eliminates  the  need  to  evaluate  the  effective  matrix  stress  rate,  Ym,  as  was 
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required  in  Equation  (6).  The  value  of  the  effective  matrix  stress  can  be  obtained  from  the  yield 
function  (Equation  (1)): 


(2+p2)cJ2  +^-ll 

sTpj 


(C.  9 ) 


The  yield  function  described  by  Equation  (1)  is  treated  as  a  plastic  potential.  Implicitly  assumed 
is  that  the  loading  and  unloading  are  controlled  by  the  BP  model  through  the  matrix  equivalent 
plastic  strain. 

Examination  of  Equation  (B.7)  shows  that  Df  (-  (D%) 2  )  is  ultimately  a  function  of 

the  same  four  integration  variables  (the  variables  input  from  EPIC  are  considered  constant  during 
the  EPIC  time  step).  The  BP  state  variable  Z  is  only  a  function  of  Wp  through  Equation  (C.7), 

and  the  term  3  J2  (=V^)  is  a  function  of  f,  P,  and  Sy  through  Equation  (C.9).  The  variables 

f,  P,  and  Sy  are  in  turn  functions  of  e?j  through  Eqs.  (24)  to  (26)  and  Eqs.  (C.2)  to  (C.6).  Since 

the  variables  P  and  Sy  are  also  direct  functions  of  fn  (through  the  degraded  modulus),  and  the 
fn  is  in  turn  a  function  of  P  and  Sy  through  Equations  (C.8)  and  (C.9),  an  efficient  convergent 
iterative  scheme  was  developed  to  solve  fn  simultaneously  with  P  and  Sy.  Since  the  temperature 
(T)  in  Equation  (B.13)  and  the  pressure  (Ps)  in  Equation  (12)  are  both  direct  functions  of  the 
internal  energy,  Equation  (38)  must  also  be  integrated. 

This  rearrangement  of  the  governing  equations  has  resulted  in  fewer  differential  equations 
which  are  relatively  less  stiff.  Thus,  the  numerical  integration  problem  was  reduced  to  the 

solution  of  the  rate  equations  for  D%,  Wp,  I,  and  The  corresponding  rate  equations  were 

solved  using  the  DIRK  scheme.  The  solution  stability,  accuracy,  and  uniqueness,  with  respect 
to  different  finite  element  mesh  sizes  and  time  steps,  were  demonstrated  through  one-dimensional 
plate  impact  solutions  in  Reference  [C.l]  and  through  two-dimensional  solutions  of  several 
application  problems. 
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Appendix  D 

Diagonally  Implicit  Runge  Kutta  Scheme 


The  diagonally  implicit  Runge-Kutta  integration  is  given  by  the  formula  [Dl], 

icjn)  =  f  (tn  +  Cjb,  yn  +  hZaijK(jn) )  ,  i=l(l)g  (D«D 

yn.i  =y„+AIbiKj”’  (d-2) 

{&)n  =  tun.y„)  (D.3) 

h  is  the  time  step  and  yn  is  the  vector  of  the  integrated  variables  at  time,  t^  Note  that  Ks  is  also 
a  function  of  itself,  which  requires  an  iterative  corrector  procedure.  That  is,  the  calculated  value 
of  Kj  is  re-substituted  into  the  nonlinear  equation  (D.l)  until  convergence  is  achieved  within  an 
error  tolerance.  If  convergence  is  not  achieved  after  a  certain  number  of  re-substitutions  of  Kj, 
then  the  time  step,  h,  is  reduced.  Convergence  of  K;  is  assured  if  a  small  enough  time  step  is 
chosen.  In  the  RDG  model  implementation,  if  convergence  is  achieved  within  a  very  small  error 
tolerance,  the  time  step  h  is  gradually  increased  to  improve  computation  time.  This  scheme  was 
implemented  into  the  shock  wave  propagation,  finite  element  code,  EPIC-2. 

The  EPIC-2  host  code  is  based  on  the  time-centered  integration  approach  with  the  time 
step  controlled  by  the  Courant  stability  criterion;  therefore,  the  time  steps  usually  remain  small. 
Thus,  the  lowest  order  DIRK,  (q=l),  will  be  adequate  for  stable  and  accurate  solutions.  For  the 
values,  q=l,  an=l/2,  bj=l,  and  Cj=l/2,  the  DIRK  method  becomes  second-order  accurate.  This 
special  case  is  equivalent  to  Bass  and  Oden’s  [D2]  recommendation  of  their  Euler  predictor- 
trapezoidal  corrector  scheme.  Saliba  [D3]  has  effectively  proved  this  special  case  as  stable  and 
accurate  in  conjunction  with  the  mechanical  equilibrium  equations  for  quasi-static  problems. 

With  the  DIRK  method  applied  to  the  ODE  of  any  viscoplastic  formulation,  a  single 
corrector  step  was  found  sufficient  for  accurate  calculations  most  of  the  time.  However,  during 
the  stiff  phase  of  the  solutions,  such  as  the  unloading,  the  DIRK  time  steps  require  significant 


D-l 


reduction  even  with  several  re-substitutions  of  Kj.  Maximum  efficiency  during  the  stiff  phases 
was  achieved  as  follows.  If  after  three  re-substitutions  of  Kj  the  relative  error  is  not  within  two 
percent,  the  DIRK  time  step  is  reduced  by  one-tenth.  The  process  continues  until  the  two  percent 
error  is  obtained.  If,  however,  after  three  re-substitutions  of  Kj  the  relative  error  is  within  two- 
tenth  of  a  percent,  the  DIRK  time  step  is  increased  in  an  arbitrary  manner  by  33%  for  the  next 
time  step.  This  process  continues  until  the  end  of  the  EPIC  time  step  is  reached. 
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Appendix  E 

Modified  TCK  Model  for  Ceramic  Materials 


The  following  pages  were  published  as  a  journal  article  entitled,  "Impact  Damage  Model 
for  Ceramic  Materials,"  J.  Appl.  Phys .,  Vol.  66,  No.  8,  15  October,  1989. 
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The  fracture  process  in  ceramic  materials  upon  impact  loading  is  complex  in  nature.  Most 
often,  the  brittle  ceramic  deforms  inelastically  due  to  microcracking  under  shock 
(compression)  loading.  At  high-velocity  impact,  the  shock  generated  microcracks  rapidly 
open  and  extend  under  subsequent  tension  (due  to  the  release  waves  from  the  stress-free 
boundaries)  leading  to  complete  pulverization  of  the  ceramic  materials.  The  main  objective  of 
this  paper  is  to  model  the  damage  process  in  ceramics  due  to  impact  loading.  A 
computationally  oriented  continuum  damage-based  constitutive  model  is  considered.  Several 
modifications  incorporating  strain  rate  and  damage  effects  on  the  compressive  strength  have 
been  introduced  into  the  model.  Results  are  presented  in  terms  of  numerical  simulations  of  a 
plate-impact  test  configuration.  Effects  of  the  model  parameters  on  the  compressive  strength 
and  spall  strength  are  described.  The  proposed  damage  model  has  been  used  successfully  to 
match  the  measured  stress  history  from  a  plate-impact  experiment  on  AD-85  (85%  aluminum 
oxide)  ceramic  target  material. 


I.  INTRODUCTION 

Ceramic  materials  usually  exhibit  very  large  compres¬ 
sive  strength  when  compared  to  metals.  The  compressive 
strength  is  found  to  increase  significantly  under  dynamic 
loading  conditions.  The  increase  in  strength  is  usually  attrib¬ 
uted  to  the  combined  effects  of  confining  pressure  and  strain 
rate.  Due  to  their  high  strength  and  enhanced  high-tempera¬ 
ture  properties,  ceramic  materials  have  often  been  employed 
as  armor  elements  and  turbine  engine  components,  respec¬ 
tively.  The  increased  use  of  ceramic  as  armor  material  neces¬ 
sitates  a  thorough  understanding  of  ceramic  material  behav¬ 
ior  under  impact  loading  conditions.  Initially,  in  armor 
design,  the  response  of  the  ceramic  armor  or  target  was  inter¬ 
polated  or  sometimes  even  extrapolated  from  empirical 
curves  constructed  from  experimental  data.  With  the  advent 
of  increasing  computer  capabilities,  the  experimental  design 
and  amnor  material  selection  are  often  guided  by  computer 
simulations  of  the  impact  situations.  General  purpose  shock- 
wave-propagation  finite-element-diffcrence  codes  have  been 
widely  used  in  the  initial  estimations  of  the  armor  perfor¬ 
mance.  The  accuracy  and  predictability  of  the  numerical  cal¬ 
culations  depend  on  the  realistic  description  of  the  ceramic 
material  through  appropriate  constitutive  and  damage  mod¬ 
els  in  the  code. 

Several  investigators1-*  have  modeled  the  behavior  of 
rock-type  brittle  solids.  In  these  models,  the  effects  of  dam¬ 
age  on  the  constitutive  behavior  were  introduced  through 
degraded  elastic  moduli.  Costin  and  Stone2  and  Horii  and 
Nemat-Nasscr’  developed  constitutive  models  for  brittle 
materials  under  static  loading  conditions  where  strain-rate 
effect  on  the  damage  process  is  not  considered.  These  models 
are  capable  of  modeling  different  failure  modes,  such  as  the 
faulting  and  splitting,  and  describe  damage-induced  aniso¬ 
tropic  material  behavior.  In  general,  these  models  are  com¬ 
putationally  demanding  and  mathematically  complex,  espe¬ 
cial!)'  the  model  by  Horii  and  Ncmat-Nasser.  The 
microphysical  model  of  Margolin*  considers  both  frictional 


effects  due  to  closed  cracks  and  strain-rate  effects.  This  mod¬ 
el  treats  damage  as  a  scalar  variable  and  is  based  on  degraded 
elastic  moduli.  Recently,  Taylor,  Chen,  and  Kuszmaul5  pre¬ 
sented  a  constitutive  or  damage  model  to  describe  the  behav¬ 
ior  of  oil  shale  under  explosive  loading  conditions.  Follow¬ 
ing  the  approach  of  Grady  and  Kipp,'  this  relatively  simple 
model  incorporates  strain-rate  effects  for  tensile  loading. 
The  inelastic  response  of  rocks  under  shock  (compression) 
loading  is  treated  simply  as  an  elastic-perfectly-plastic  re¬ 
sponse.  The  compressive  strength  K is  treated  as  a  constant 
independent  of  damage,  pressure,  and  strain  rate.  Rajen¬ 
dran,  Kroupa,  and  Brar'*  implemented  the  model  by  Taylor, 
Chen,  and  Kuszmaul  (TCK  model)  in  a  two-dimensional 
finite-element  code,  EPIC-2,  and  simulated  a  plate-impact  ex¬ 
perimental  configuration  using  this  model.  Recently,  Rajen¬ 
dran  and  Cook,,0in  a  review  report,  cited  several  experimen¬ 
tal  results  that  show  the  effects  of  pressure  and  strain  rate  on 
the  compressive  strength  of  various  ceramics. 

Hence,  several  modifications  to  the  original  model  are 
required  for  a  realistic  description  of  ceramic-type  brittle 
materials.  In  the  present  paper,  the  application  of  the  modi¬ 
fied  TCK  model  to  describe  ceramic  material  behavior  un¬ 
der  impact  loading  conditions  is  considered.  In  the  next  sec¬ 
tion,  the  original  TCK  model  is  described.  A  strain-rate  and 
damage-dependent  compressive  strength  model  is  proposed 
and  incorporated  into  the  constitutive  formulation.  For  ten¬ 
sile  loading,  the  original  TCK  model  is  retained,  except  that 
the  unloading  and  reloading  paths  are  different  in  the  pres¬ 
ent  approach.  Since  the  proposed  model  for  ceramic  materi¬ 
als  is  useful  for  computer  code  applications,  we  used  an  im¬ 
proved  and  efficient  numerical  algorithm  for  solving  the 
governing  constitutive  equations.  A  second-order  diagonal¬ 
ly  implicit  Runge-Kutta  (DIRK)  method  is  employed  in 
the  numerical  solution.  The  effects  of  time  step  and  grid  sen¬ 
sitivity  on  the  stability  and  accuracy  of  the  solutions  are  veri¬ 
fied.  Using  the  proposed  constitutive  or  damage  model  the 
response  of  AD-85  (85%  aluminum  oxide)  ceramic  under 
plate-impact  loading  (uniaxial  strain)  conditions  was  mod- 
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clod.  The  experimentally  determined  stress  history  matched 
extremely  well  with  the  numerical  simulation. 


II.  TAYLOR-CHEN-KUSZMAUL  MODEL 

Under  dynamic  compressive  loading  conditions,  rock- 
and  concretelike  materials  deform  inelastically  similar  to  ce¬ 
ramic  materials.  The  inelastic  response  of  these  brittle  mate¬ 
rials  is  usually  attributed  to  the  nucleation,  growth,  and  co¬ 
alescence  of  microcracks.  The  original  TCK  model5  treats 
the  inelastic  behavior  under  compression  as  an  elastic-per- 
fectly-plastic  response.  This  elastic-perfectly-plastic  as¬ 
sumption  simplifies  the  model  formulation  for  compressive 
loading,  and  therefore  a  Prandtl-Reuss  or  Prandtl-Mises 
incremental-plasticity-theory-based  model  is  adequate  to 
describe  the  inelastic  behavior. 

A  crack-growth-based  damage  model  has  been  used  to 
describe  the  inelastic  behavior  under  tension.  The  damage- 
model  formulation  considers  the  growth  and  interaction  of 
randomly  distributed  penny-shaped  microcracks.  The  dam¬ 
age  nucleation  is  treated  by  a  Weibull  distribution.  The  num¬ 
ber  of  cracks  that  nucleate  upon  loading  is  given  by 

N=kC,  (1) 

where  A  is  the  number  of  flaws  per  unit  volume,  A:  and  m  are 
model  parameters,  and  eu  is  the  volumetric  strain.  The  above 
equation  can  be  written  in  terms  of  pressure  by  replacing  e„ 
by  P  /3 K,  where  P  is  the  pressure,  and  AT  is  the  bulk  mod  ulus. 

The  growth  and  coalescence  of  microcracks  degrade  the 
bulk  and  shear  moduli  or  the  Poisson’s  ratio  ( v) .  To  describe 
the  effects  of  crack  volume  (or  damage)  on  these  elastic 
moduli,  the  derivation  of  Budiansky  and  O’Connell"  was 
employed.  The  corresponding  expressions  are 


(2) 


v  =  v(  1  —  ^Cd). 


(3) 


The  above  expression  for  v  is  a  simplified  equivalent  form  of 
Budiansky  and  O’Connell.  K  and  vare  the  degraded  elastic 
properties  and  Cd  is  a  crack  density  parameter.  Equations 
(2)  and  (3)  are  not  sufficient  to  determine  the  three  un¬ 
knowns,  v,  K,  and  Cd.  We  need  an  additional  equation  and 
this  can  be  obtained  from  the  definition  of  crack  density,  Cd . 

The  crack  density  Cd  represents  the  volume  fraction  of 
penny-shaped  voids.  The  average  crack  density  is  defined  by 
the product  of the  number  of  cracks  per  unit  volume  (A)  and 
the  average  volume  of  the  cracks  ( ~a}),  and  is  given  by 

Cd=m\  (4) 

where  "0"  is  a  crack  geometry  parameter. 

Taylor  and  co- workers5  assumed  that  the  average  crack 
size  0  is  proportional  to  the  size  of  the  average  fragments, 
and  utilized  the  expression  derived  by  Kipp  and  Grady"  for 
fragmentation.  The  corresponding  expression  for  the  crack 
size  is 


1  /  yj20Klc 

2 


r 


(5) 


where  p  is  the  density  of  the  material,  K,c  is  the  static  frac¬ 
ture  toughness,  c  is  the  sound  speed,  and  the  fm,A  is  the 


applied  strain  rate  at  fracture.  By  combining  Eqs.  ( 1  )-(5), 
the  expression  for  Cd  is 


Q  = 


5  k 
2  (3AO 


(6) 


Model  formulation  will  be  completed  once  an  expres¬ 
sion  for  damage  is  obtained.  The  following  relationship  be¬ 
tween  damage  and  the  bulk  modulus  is  assumed: 


K  =  K{\  —  D). 


(7) 


The  bulk  modulus  is  assumed  to  degrade  linearly  with 
increasing  damage.  Also,  by  comparing  Eqs.  (2)  and  (7), 
we  can  describe  damage  in  terms  of  the  crack  density  as 


l-v>\ 

9  \ 

1  —  IvJ 

The  original  bulk  modulus  K  is  reduced  with  respect  to  dam¬ 
age  in  a  linear  manner  as  shown  in  Eq.  (7).  When  there  is  no 
damage  ( D  =  0)  the  bulk  modulus  remains  the  same  and 
K  =  K.  Upon  total  damage  (Z)  =  1)  the  material  can  no 
longer  carry  any  tensile  loading  due  to  the  complete  loss  of 
stiffness.  The  constitutive  equation  for  the  brittle  material 
can  now  be  written  as  sum  of  the  bulk  and  deviatoric  parts  of 
the  stress  components, 

au  =  m  1  -  D)ekkS0  +  2G(  1  -  D)eiJt  (9) 

where  e,t  are  the  deviatoric  part  of  the  strain  components. 
When  D  =  0,  the  Hook’s  law  for  linear  elastic  behavior  is 
recovered  from  Eq.  ( 9 ) .  The  rate  equations  for  D  and  Cd  can 
be  derived  from  Eqs.  (6)  and  (7).  The  six  components  of  the 
stress  tensor,  damage  parameter  D,  and  the  density  param¬ 
eter  Cd  are  the  unknowns  in  the  model  formulation.  These 
eight  unknowns  can  be  determined  by  solving  the  eight  rate 
equations  incrementally. 


III.  NUMERICAL  SIMULATIONS 

The  TCK  model  was  implemented  into  the  EPIC-2  finite- 
element  code."  Initially,  we  developed  a  fourth-order 
Runge-Kutta  scheme  to  solve  the  governing  equations  in¬ 
crementally.  However,  an  improved  scheme  based  on  a  sec¬ 
ond-order  implicit  Runge-Kutta  method  was  developed  lat¬ 
er  for  numerical  stability  and  accuracy.  For  model 
validation,  a  plate-impact  test  on  a  ceramic  target  was  con¬ 
sidered. 


A.  Plate-impact  test  description 

A  schematic  representation  of  the  plate-impact  test  is 
shown  in  Fig.  1 .  A  copper  flyer  disk  of  50  mm  diameter  and 
2.5  mm  thickness  impacts  a  AD-85  ceramic  plate  of  thick- 
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FIG.  I.  A  schematic  representation  of  the  plate-impact  experiment. 
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ness  S.S  mm.  The  ceramic  plate  was  backed  by  a  12-mm- 
thick  PMMA  (polymelhy)  methacrylate).  The  stress  his¬ 
tory  at  the  interface  of  AD-85  and  PMMA  was  recorded 
using  a  Manganin  gauge.  For  an  impact  velocity  of  570  m/s, 
the  stress-time  profile  from  the  Manganin  gauge  is  shown  in 
Fig.  2.  The  interpretation  of  this  stress  profile  is  as  follows: 
The  point  "A  "  on  the  profile  corresponds  to  the  elastic 
strength  of  AD-85  under  one-dimensional  strain  condition 
and  it  is  known  as  Hugoniot  elastic  limit  ( HEL)  of  the  mate¬ 
rial.  The  profile  between  A  and  the  peak  stress  indicated  by 
point  B  represents  the  inelastic  (some  call  it  “plastic")  part 
of  the  stress  w  ave.  Point  C  corresponds  to  the  arrival  of  the 
release  wave  which  originated  at  the  rear  surface  of  the  flyer. 
The  profile  between  points  C  and  D  is  due  to  elastic  unload¬ 
ing  of  the  stress  release.  The  release  characteristics  are  quite 
complex  to  interpret.  Unlike  in  metals,  most  ceramics  exhib¬ 
it  a  lack  of  a  deterministic  spall  signal,  and  therefore  some 
kind  of  qualitative  interpretations  only  can  be  made.  The 
slow  stress  release  in  steps  between  points  D  and  E  indicates 
zero  spall  strength  of  the  material.  The  interpretations  of 
stress  gauge  profiles  from  plate-impact  tests  on  ceramic  ma¬ 
terials  were  reviewed  in  detail  by  Bless.14 

B.  Effect  of  time  step 

In  the  next  stage,  we  simulated  the  plate-impact  test 
configuration  using  the  epic-2  code  and  the  TCK  model.  A 
one-dimensional  strain  condition  is  assumed  in  the  simula¬ 
tion.  The  stress  history  at  the  interface  between  the  ceramic 
target  and  the  PMMA  (backup)  plate  was  obtained  from 
the  numerical  simulation.  Several  computer  runs  were  made 
with  different  time  steps  in  each  run.  The  stress  histories 
from  these  runs  are  compared  in  Fig.  3.  It  can  be  seen  that 
the  time  step  had  no  influence  on  the  stress  between  points  A 
and  E in  the  figure.  This  is  because  the  constitutive  model  for 
compression  follows  a  simple  elastic-plastic  approach  with¬ 
out  any  damage  calculation  and  no  unloading  is  involved. 
However,  beyond  point  E,  the  stress  calculations  involve 


FIG.  2.  Stress  history  in  PMMA  from  the  plate-impact  test  on  AD-85. 
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FIG.  3.  The  effect  of  time  step  on  (he  numerical  solutions. 


damage  nucleation,  growth,  and  degradation  of  stiffness  and 
strength  leading  to  target  spallation.  The  numerical  solution 
exhibited  both  grid  and  time-step  dependency,  indicating 
nonunique  and  unstable  solution,  during  the  stress  release 
(due  to  spallation).  To  investigate  this  aspect  further,  the 
unloading  and  reloading  nature  of  the  model  formulation 
was  critically  evaluated.  We  simulated  the  hypothetical 
loading,  unloading  and  reloading  conditions  shown  in  Fig.  4. 
In  the  TCK  model,  unloading  is  along  the  undamaged  mod¬ 
ulus  (between  points  A  and  B)  while  reloading  (between 
points  B  and  C)  is  along  the  damaged  modulus  as  show  n  in 
the  figure.  When  reloading  occurs  at  different  points  (B  ', 
B  ’ )  the  pressure  reaches  different  levels  for  the  same  strain, 
indicating  nonunique  solutions.  Brittle  solids  usually  do  not 
exhibit  any  gross  plastic  deformation,  and  therefore  strains 


0000  0  002  0004  0006  0.008  0  010 


STRAIN 

FIG.  4.  Nonunique  solutions  due  to  unloading  along  the  undamaged  modu¬ 
lus  and  reloading  along  the  damaged  modulus. 
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arc  not  accumulated  upon  unloading.  The  unloading  occurs 
along  the  degraded  modulus  toward  the  origin  of  the  uniax¬ 
ial  stress-strain  curve.  Therefore,  in  the  model,  we  imple¬ 
mented  unloading  and  reloading  along  the  current  damaged 
modulus  until  damage  reinitiated  for  further  analysis.  Addi¬ 
tional  plate-impact  test  simulations  were  conducted  with 
this  feature.  The  solutions  obtained  using  different  time 
steps  yielded  almost  the  same  stress  history  as  can  be  seen 
from  Fig.  5. 


C.  Effects  of  model  parameters 

There  are  four  model  parameters  (without  including 
the  bulk  and  shear  moduli  of  the  intact  ceramic ) :  the  nuclea- 
tion  parameters  k  and  m,  the  fracture  toughness  K,c,  and  the 
compressive  yield  strength  Y.  The  compressive  yield 
strength  can  be  directly  obtained  from  the  stress  (oHEL  )  at 
HEL.  The  corresponding  relationship  between  Y  and  <7HEL 
is  given  by 


^HEL 

K/lG+\ 


(10) 


The  values  of  k  and  m  have  to  be  calibrated  based  on  the 
model’s  ability  to  reproduce  the  experimental  stress  history. 
A  value  for  K,c  can  be  determined  from  a  standard  fracture 
toughness  test.  However,  for  better  understanding  of  these 
three  parameters  and  their  effects  on  the  spall  signal  (be¬ 
tween  points  D  and  E  of  Fig.  2 )  of  the  stress  history,  several 
numerical  simulations  were  carried  out  by  varying  one  pa¬ 
rameter  at  a  time. 

In  Fig.  6,  the  effects  of  the  damage  nucleation  parameter 
“k  ”  are  shown  for  m  =  6  and  K,c  —  3.0  (MPa)  >[m.  The 
Stress  histories  at  the  stress  gauge  location  for  Ar  =  10’°,  10”, 
and  10'4  are  plotted.  For  k  <  lO’’",  the  damage  (micro¬ 
cracks)  level  is  relatively  low  and  in  the  simulation  the  ce¬ 
ramic  target  did  not  fracture.  The  stress  history  follows  the 


FIG.  5.  Effect  of  time  step  on  the  stable  numerical  solutions  due  to  unload¬ 
ing  and  reloading  along  the  damaged  modulus. 


FIG.  6.  Effect  of  damage  nucleation  parameter  k  on  the  spall  signal. 


solid  line  AB,  which  indicates  complete  unloading  of  the 
stress  due  to  the  release  wave  from  the  back  surface  of  the 
flyer  plate.  When  the  value  of  k  is  increased  to  10”  and  10'4, 
the  ceramic  target  fractured  by  spalling  and  the  slow  step 
release  of  stress  between  A  'B '  or  A  "B  "  clearly  shows  the  sa¬ 
lient  features  of  an  experimentally  obtained  stress  history 
(see  Fig.  2).  A  higher  value  of  k  represents  a  larger  number 
of  cracks  nucleating,  and  therefore  earlier  arrival  (at  point 
A  ’ )  of  the  spall  signal. 

Figure  7  demonstrates  how  m  also  controls  the  damage 
nucleation  level.  In  these  simulations,  a  value  of  10”  for  k 
and  3.0  (MPa )Jm  for  K,c  are  employed.  Stress  histories 
were  obtained  for  m  -  4,  5,  6,  and  7.  When  the  value  of  m 
was  higher,  the  damage  nucleation  was  lower  [note  that  in 
Eq.  (1),  e„<l),  and  therefore  the  spall  signal  arrival  oc¬ 
curred  at  a  later  time.  For  values  of  m  less  than  5,  the  spall 
signal  was  relatively  unaffected,  indicating  complete  pulver¬ 
ization  of  the  target.  Unlike  in  metals  where  the  spall  is  limit- 


FIG.  7.  Effect  of  the  damage  nucleation  parameter  m  on  the  spall  signal. 
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ed  to  a  single  plane,  in  ceramics  spall  occurs  over  a  zone  as 
can  be  clearly  seen  from  this  figure.  This  is  confirmed  by  the 
experimental  observation  of  spall  zones  in  polycrystalline 
ceramics  by  Bless,  Yaziv,  and  Rosenberg. 15  Most  of  the  tar¬ 
get  experienced  a  critical  damage  level  indicating  pulveriza¬ 
tion.  It  is  clear  from  these  simulations  (effects  of  k  and  m  on 
the  spall  signal )  that  the  spall  arrival  time  in  the  experiment 
can  be  used  to  estimate  the  values  of  k  and  m.  An  arbitrary 
value  of  6  for  m  seems  to  work  well  in  the  experimental 
calibration.  The  value  of  k  can  be  adjusted  to  match  the  spall 
arrival  time. 

As  a  final  exercise,  the  value  of  K,c  was  varied  (m  =  6 
and  k  =  1027)  and  stress  history  at  the  gauge  location  was 
calculated.  The  results  indicated  less-pronounced  effects  of 
the  fracture  toughness  on  the  spall  processes  as  can  be  seen 
from  Fig.  8.  Shock  loading  induces  microcracks  in  most  ce¬ 
ramics.  These  microcracks  often  nucleate  from  preexisting 
flaws,  thus  degrading  subsequent  tensile  properties.  There¬ 
fore,  a  small  tensile  load  is  sufficient  to  propagate  the  cracks. 
Hence  K,c  may  not  have  significant  effect  on  the  spall  re¬ 
lease. 

D.  Modeling  of  AO-85  ceramic 

We  modeled  the  plate-impact  experiment  described  ear¬ 
lier  in  this  section.  The  target  was  AD-85  ceramic.  The  fol¬ 
lowing  values  of  the  model  parameters  were  used:  K,c 
=  3  0(MPa )Jm,  k  =  1027,  m  =  6,  Y  =  4000  MPa.  Note 
that  the  value  of  K,c  is  obtained  from  open  literature  as  the 
fracture  toughness  value  for  AD-85. 

A  comparison  between  the  stress  history  calculated  in 
the  simulation  and  the  experiment  is  shown  in  Fig.  9.  It  can 
be  seen  that  the  stress  profile  between  points  4  (at  HEL)  and 
B  ( the  peak  stress)  does  not  match  well.  The  ramping  nature 
of  the  stress  versus  time  curve  between  these  two  points  in 
the  experiment  may  be  speculated  as  due  to  some  sort  of 
microdamaging  mechanism  which  is  a  rate-controlled  pro¬ 
cess.  In  metals,  this  ramping  is  usually  attributed  to  the 


FIG.  8.  Effect  of  K,c  on  the  spall  signal. 
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FIG.  9.  Calculated  stress  history  using  a  strain-rate-independent  model  for 
compression. 


visco-plasticity  (strain-rate  sensitivity  of  the  flow  stress). 
The  behavior  of  ceramic  under  compression  is  assumed  to  be 
elastic-perfectly-plastic  in  the  numerical  simulation  using 
the  original  TCK.  model.  Since  the  strain-rate  dependency 
was  absent,  the  simulated  stress  history  showed  the  typical 
‘‘knee’’  at  the  HEL.  The  “plastic"  wave  profile  was  without 
any  ramping  which  is  typically  observed  in  strain-rate-inde¬ 
pendent  metals. 

For  improved  modeling  of  the  behavior  of  ceramic  un¬ 
der  compression,  the  following  form  for  the  compressive 
strength  is  considered l0: 

Y=  r,(l  -Finn e)(l-Z»,  (II) 

where  Y,  is  the  static  strength  and  B  is  the  strain-rate  sensi¬ 
tivity  parameter,  e  is  the  equivalent  plastic  strain  rate  and  D 
is  the  damage.  The  damage  growth  is  assumed  to  be  zero 
under  compression.  However,  the  compression  strength  is 
reduced  after  accumulating  damage  under  tension.  We  in¬ 
corporated  this  compression  model  into  the  formulation  and 
resimulated  the  plate-impact  experiment.  The  value  of  Y,  is 
determined  from  static  compressive  strength  as  2000  MPa. 
By  knowing  the  T  that  corresponds  to  the  HEL,  a  value  of 
0.175  was  obtained  for  the  parameter  B.  The  calculated 
stress  history  at  the  gauge  location,  based  on  the  new  model, 
compared  extremely  well  with  the  experimental  data  as  can 
be  seen  from  Fig.  10. 

IV.  SUMMARY 

In  general,  the  behavior  of  ceramic  under  impact  load¬ 
ing  conditions  is  complex.  Since  a  recovery  test  on  ceramic  is 
difficult  to  perform  especially  at  high-velocity  impact 
(above  the  HEL),  many  of  the  deformation  and  damage 
mechanisms  are  speculative  in  nature.  The  free-surface  ve¬ 
locity  history  of  a  ceramic  target  and  the  stress  history  at  the 
interface  of  ceramic  plate  and  the  backup  plate  are  the  only 
measured  quantities  which  we  rely  upon  in  the  ceramic  ma¬ 
terial  modeling.  At  present,  the  capabilities  of  any  advanced 
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FIG.  10.  Comparison  belween  the  experimental  and  the  modified 
TCK  model. 

engineering  or  physically  based  model  are  validated  only 
through  the  model’s  ability  to  match  the  various  salient  fea¬ 
tures  of  the  experimentally  measured  data. 

A  constitutive  or  damage  model  for  ceramic  material 
under  impact  loading  conditions  was  presented.  This  model 
employed  strain-rate-dependent  equations  for  describing  the 
inelastic  behavior  of  ceramic  under  both  compression  and 
tension.  The  tension  model  was  originally  developed  for  de¬ 
scribing  fragmentation  of  rocks  and  oil  shale.  A  more  realis¬ 
tic  description  of  unloading  and  loading  paths  was  incorpo¬ 
rated  in  the  model  for  unique  and  stable  solutions.  The 
model  was  then  implemented  into  a  shock-wave  propaga¬ 
tion,  finite-element  code.  The  stability  of  the  numerical  solu¬ 
tion  was  demonstrated.  A  model  parameter  evaluation 
scheme  is  presented  in  terms  of  the  sensitivity  study.  We 
have  successfully  matched  the  experimental  stress  history 
obtained  from  a  plate-impact  test  on  a  AD-85  ceramic.  The 


impact  behaviors  of  different  ceramics  can  be,  in  general, 
very  different.  Depending  on  the  impact  velocity  and  ceram¬ 
ic  constituents,  deformation  and  fracture  mechanisms  can 
vary  significantly  from  one  ceramic  to  another.  The  capabili¬ 
ties  of  the  proposed  model  to  describe  other  ceramic  materi¬ 
als  are  yet  to  be  validated. 
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